MODULE COMP_FUNCTIONS

! I/O + OS functions

USE PRECISION_PARAMETERS 
IMPLICIT NONE 
CHARACTER(255), PARAMETER :: funcid='$Id: func.f90 6239 2010-05-17 15:14:52Z mcgratta $'
CHARACTER(255), PARAMETER :: funcrev='$Revision: 6239 $'
CHARACTER(255), PARAMETER :: funcdate='$Date: 2010-05-17 16:14:52 +0100 (Mon, 17 May 2010) $'
 
CONTAINS

SUBROUTINE GET_REV_func(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') funcrev(INDEX(funcrev,':')+1:LEN_TRIM(funcrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') funcdate

END SUBROUTINE GET_REV_func

REAL(EB) FUNCTION SECOND()  ! Returns the CPU time in seconds.
!$ USE OMP_LIB
REAL(FB) CPUTIME
CALL CPU_TIME(CPUTIME)
SECOND = CPUTIME
!$ SECOND = OMP_GET_WTIME()
END FUNCTION SECOND


REAL(EB) FUNCTION WALL_CLOCK_TIME()  ! Returns the number of seconds since January 1, 2000, including leap years

! Thirty days hath September,
! April, June, and November.
! February has twenty-eight alone;
! All the rest have thirty-one,
! Excepting Leap-Year, that's the time
! When February's days are twenty-nine.

INTEGER :: DATE_TIME(8),WALL_CLOCK_SECONDS
CHARACTER(10) :: BIG_BEN(3)
! X_1 = common year, X_2 = leap year
INTEGER, PARAMETER :: S_PER_YEAR_1=31536000,S_PER_YEAR_2=31622400,S_PER_DAY=86400,S_PER_HOUR=3600,S_PER_MIN=60
INTEGER, PARAMETER, DIMENSION(12) :: ACCUMULATED_DAYS_1=(/0,31,59,90,120,151,181,212,243,273,304,334/), & 
                                     ACCUMULATED_DAYS_2=(/0,31,60,91,121,152,182,213,244,274,305,335/)
INTEGER :: YEAR_COUNT

CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
WALL_CLOCK_SECONDS = 0._EB
DO YEAR_COUNT=2001,DATE_TIME(1)
   !Leap year if divisible by 4 but not 100 unless by 400 (1900 no, 1904  yes, 2000 yes)
   IF (MOD(YEAR_COUNT,4)==0 .AND. (MOD(YEAR_COUNT,100)/=0 .OR. MOD(YEAR_COUNT,400)==0)) THEN
      WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_YEAR_2
   ELSE
      WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_YEAR_1
   ENDIF
ENDDO 
IF (MOD(DATE_TIME(1),4)==0 .AND. (MOD(DATE_TIME(1),100)/=0 .OR. MOD(DATE_TIME(1),400)==0 )) THEN
   WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_DAY*(ACCUMULATED_DAYS_2(DATE_TIME(2))+DATE_TIME(3))
ELSE
   WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_DAY*(ACCUMULATED_DAYS_1(DATE_TIME(2))+DATE_TIME(3))
ENDIF
WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS +  S_PER_HOUR*DATE_TIME(5) + S_PER_MIN*DATE_TIME(6) + DATE_TIME(7)
WALL_CLOCK_TIME    = WALL_CLOCK_SECONDS + DATE_TIME(8)*0.001_EB

END FUNCTION WALL_CLOCK_TIME


SUBROUTINE SHUTDOWN(MESSAGE)  

! Stops the code gracefully after writing a message

USE GLOBAL_CONSTANTS, ONLY: LU_ERR,FN_OUTPUT,LU_OUTPUT,CHID
CHARACTER(*) MESSAGE
LOGICAL :: EX

WRITE(LU_ERR,'(/A)') TRIM(MESSAGE)

! For those running FDS on a PC, ensure that error message gets written to the .out file, too.

INQUIRE(FILE=FN_OUTPUT,EXIST=EX)
IF (.NOT.EX) OPEN(LU_OUTPUT,FILE=TRIM(CHID)//'.out',STATUS='REPLACE',FORM='FORMATTED')

WRITE(LU_OUTPUT,'(/A)') TRIM(MESSAGE)

STOP

END SUBROUTINE SHUTDOWN


SUBROUTINE GET_INPUT_FILE ! Read the argument after the command
USE GLOBAL_CONSTANTS, ONLY: FN_INPUT
IF (FN_INPUT=='null') CALL GETARG(1,FN_INPUT)
END SUBROUTINE GET_INPUT_FILE


INTEGER FUNCTION GET_FILE_NUMBER()
USE GLOBAL_CONSTANTS, ONLY: MYID,FILE_COUNTER
FILE_COUNTER(MYID) = FILE_COUNTER(MYID) + 1
GET_FILE_NUMBER = FILE_COUNTER(MYID)
END FUNCTION GET_FILE_NUMBER



SUBROUTINE CHECKREAD(NAME,LU,IOS)
 
INTEGER :: II,IOS
INTEGER, INTENT(IN) :: LU
CHARACTER(4), INTENT(IN) :: NAME
CHARACTER(80) TEXT
 
IOS = 1
READLOOP: DO
   READ(LU,'(A)',END=10) TEXT
   TLOOP: DO II=1,72
      IF (TEXT(II:II)/='&' .AND. TEXT(II:II)/=' ') EXIT TLOOP
      IF (TEXT(II:II)=='&') THEN
         IF (TEXT(II+1:II+4)==NAME) THEN
            BACKSPACE(LU)
            IOS = 0
            EXIT READLOOP
         ELSE
            CYCLE READLOOP
         ENDIF
      ENDIF
   ENDDO TLOOP
ENDDO READLOOP
 
10 RETURN
END SUBROUTINE CHECKREAD


SUBROUTINE CHECK_XB(XB)
REAL(EB) :: DUMMY,XB(6)
INTEGER  :: I
DO I=1,5,2
   IF (XB(I)>XB(I+1)) THEN
      DUMMY   = XB(I)
      XB(I)   = XB(I+1)
      XB(I+1) = DUMMY
   ENDIF
ENDDO
END SUBROUTINE CHECK_XB


SUBROUTINE COLOR2RGB(RGB,COLOR)

! Translate character string of a color name to RGB value

INTEGER :: RGB(3)
CHARACTER(25) :: COLOR
CHARACTER(100) :: MESSAGE

SELECT CASE(COLOR)
CASE ('ALICE BLUE');RGB = (/240,248,255/)
CASE ('ANTIQUE WHITE');RGB = (/250,235,215/)
CASE ('ANTIQUE WHITE 1');RGB = (/255,239,219/)
CASE ('ANTIQUE WHITE 2');RGB = (/238,223,204/)
CASE ('ANTIQUE WHITE 3');RGB = (/205,192,176/)
CASE ('ANTIQUE WHITE 4');RGB = (/139,131,120/)
CASE ('AQUAMARINE');RGB = (/127,255,212/)
CASE ('AQUAMARINE 1');RGB = (/118,238,198/)
CASE ('AQUAMARINE 2');RGB = (/102,205,170/)
CASE ('AQUAMARINE 3');RGB = (/69,139,116/)
CASE ('AZURE');RGB = (/240,255,255/)
CASE ('AZURE 1');RGB = (/224,238,238/)
CASE ('AZURE 2');RGB = (/193,205,205/)
CASE ('AZURE 3');RGB = (/131,139,139/)
CASE ('BANANA');RGB = (/227,207,87/)
CASE ('BEIGE');RGB = (/245,245,220/)
CASE ('BISQUE');RGB = (/255,228,196/)
CASE ('BISQUE 1');RGB = (/238,213,183/)
CASE ('BISQUE 2');RGB = (/205,183,158/)
CASE ('BISQUE 3');RGB = (/139,125,107/)
CASE ('BLACK');RGB = (/0,0,0/)
CASE ('BLANCHED ALMOND');RGB = (/255,235,205/)
CASE ('BLUE');RGB = (/0,0,255/)
CASE ('BLUE 2');RGB = (/0,0,238/)
CASE ('BLUE 3');RGB = (/0,0,205/)
CASE ('BLUE 4');RGB = (/0,0,139/)
CASE ('BLUE VIOLET');RGB = (/138,43,226/)
CASE ('BRICK');RGB = (/156,102,31/)
CASE ('BROWN');RGB = (/165,42,42/)
CASE ('BROWN 1');RGB = (/255,64,64/)
CASE ('BROWN 2');RGB = (/238,59,59/)
CASE ('BROWN 3');RGB = (/205,51,51/)
CASE ('BROWN 4');RGB = (/139,35,35/)
CASE ('BURLY WOOD');RGB = (/222,184,135/)
CASE ('BURLY WOOD 1');RGB = (/255,211,155/)
CASE ('BURLY WOOD 2');RGB = (/238,197,145/)
CASE ('BURLY WOOD 3');RGB = (/205,170,125/)
CASE ('BURLY WOOD 4');RGB = (/139,115,85/)
CASE ('BURNT SIENNA');RGB = (/138,54,15/)
CASE ('BURNT UMBER');RGB = (/138,51,36/)
CASE ('CADET BLUE');RGB = (/95,158,160/)
CASE ('CADET BLUE 1');RGB = (/152,245,255/)
CASE ('CADET BLUE 2');RGB = (/142,229,238/)
CASE ('CADET BLUE 3');RGB = (/122,197,205/)
CASE ('CADET BLUE 4');RGB = (/83,134,139/)
CASE ('CADMIUM ORANGE');RGB = (/255,97,3/)
CASE ('CADMIUM YELLOW');RGB = (/255,153,18/)
CASE ('CARROT');RGB = (/237,145,33/)
CASE ('CHARTREUSE');RGB = (/127,255,0/)
CASE ('CHARTREUSE 1');RGB = (/118,238,0/)
CASE ('CHARTREUSE 2');RGB = (/102,205,0/)
CASE ('CHARTREUSE 3');RGB = (/69,139,0/)
CASE ('CHOCOLATE');RGB = (/210,105,30/)
CASE ('CHOCOLATE 1');RGB = (/255,127,36/)
CASE ('CHOCOLATE 2');RGB = (/238,118,33/)
CASE ('CHOCOLATE 3');RGB = (/205,102,29/)
CASE ('CHOCOLATE 4');RGB = (/139,69,19/)
CASE ('COBALT');RGB = (/61,89,171/)
CASE ('COBALT GREEN');RGB = (/61,145,64/)
CASE ('COLD GREY');RGB = (/128,138,135/)
CASE ('CORAL');RGB = (/255,127,80/)
CASE ('CORAL 1');RGB = (/255,114,86/)
CASE ('CORAL 2');RGB = (/238,106,80/)
CASE ('CORAL 3');RGB = (/205,91,69/)
CASE ('CORAL 4');RGB = (/139,62,47/)
CASE ('CORNFLOWER BLUE');RGB = (/100,149,237/)
CASE ('CORNSILK');RGB = (/255,248,220/)
CASE ('CORNSILK 1');RGB = (/238,232,205/)
CASE ('CORNSILK 2');RGB = (/205,200,177/)
CASE ('CORNSILK 3');RGB = (/139,136,120/)
CASE ('CRIMSON');RGB = (/220,20,60/)
CASE ('CYAN');RGB = (/0,255,255/)
CASE ('CYAN 2');RGB = (/0,238,238/)
CASE ('CYAN 3');RGB = (/0,205,205/)
CASE ('CYAN 4');RGB = (/0,139,139/)
CASE ('DARK GOLDENROD');RGB = (/184,134,11/)
CASE ('DARK GOLDENROD 1');RGB = (/255,185,15/)
CASE ('DARK GOLDENROD 2');RGB = (/238,173,14/)
CASE ('DARK GOLDENROD 3');RGB = (/205,149,12/)
CASE ('DARK GOLDENROD 4');RGB = (/139,101,8/)
CASE ('DARK GRAY');RGB = (/169,169,169/)
CASE ('DARK GREEN');RGB = (/0,100,0/)
CASE ('DARK KHAKI');RGB = (/189,183,107/)
CASE ('DARK OLIVE GREEN');RGB = (/85,107,47/)
CASE ('DARK OLIVE GREEN 1');RGB = (/202,255,112/)
CASE ('DARK OLIVE GREEN 2');RGB = (/188,238,104/)
CASE ('DARK OLIVE GREEN 3');RGB = (/162,205,90/)
CASE ('DARK OLIVE GREEN 4');RGB = (/110,139,61/)
CASE ('DARK ORANGE');RGB = (/255,140,0/)
CASE ('DARK ORANGE 1');RGB = (/255,127,0/)
CASE ('DARK ORANGE 2');RGB = (/238,118,0/)
CASE ('DARK ORANGE 3');RGB = (/205,102,0/)
CASE ('DARK ORANGE 4');RGB = (/139,69,0/)
CASE ('DARK ORCHID');RGB = (/153,50,204/)
CASE ('DARK ORCHID 1');RGB = (/191,62,255/)
CASE ('DARK ORCHID 2');RGB = (/178,58,238/)
CASE ('DARK ORCHID 3');RGB = (/154,50,205/)
CASE ('DARK ORCHID 4');RGB = (/104,34,139/)
CASE ('DARK SALMON');RGB = (/233,150,122/)
CASE ('DARK SEA GREEN');RGB = (/143,188,143/)
CASE ('DARK SEA GREEN 1');RGB = (/193,255,193/)
CASE ('DARK SEA GREEN 2');RGB = (/180,238,180/)
CASE ('DARK SEA GREEN 3');RGB = (/155,205,155/)
CASE ('DARK SEA GREEN 4');RGB = (/105,139,105/)
CASE ('DARK SLATE BLUE');RGB = (/72,61,139/)
CASE ('DARK SLATE GRAY');RGB = (/47,79,79/)
CASE ('DARK SLATE GRAY 1');RGB = (/151,255,255/)
CASE ('DARK SLATE GRAY 2');RGB = (/141,238,238/)
CASE ('DARK SLATE GRAY 3');RGB = (/121,205,205/)
CASE ('DARK SLATE GRAY 4');RGB = (/82,139,139/)
CASE ('DARK TURQUOISE');RGB = (/0,206,209/)
CASE ('DARK VIOLET');RGB = (/148,0,211/)
CASE ('DEEP PINK');RGB = (/255,20,147/)
CASE ('DEEP PINK 1');RGB = (/238,18,137/)
CASE ('DEEP PINK 2');RGB = (/205,16,118/)
CASE ('DEEP PINK 3');RGB = (/139,10,80/)
CASE ('DEEP SKYBLUE');RGB = (/0,191,255/)
CASE ('DEEP SKYBLUE 1');RGB = (/0,178,238/)
CASE ('DEEP SKYBLUE 2');RGB = (/0,154,205/)
CASE ('DEEP SKYBLUE 3');RGB = (/0,104,139/)
CASE ('DIM GRAY');RGB = (/105,105,105/)
CASE ('DODGERBLUE');RGB = (/30,144,255/)
CASE ('DODGERBLUE 1');RGB = (/28,134,238/)
CASE ('DODGERBLUE 2');RGB = (/24,116,205/)
CASE ('DODGERBLUE 3');RGB = (/16,78,139/)
CASE ('EGGSHELL');RGB = (/252,230,201/)
CASE ('EMERALD GREEN');RGB = (/0,201,87/)
CASE ('FIREBRICK');RGB = (/178,34,34/)
CASE ('FIREBRICK 1');RGB = (/255,48,48/)
CASE ('FIREBRICK 2');RGB = (/238,44,44/)
CASE ('FIREBRICK 3');RGB = (/205,38,38/)
CASE ('FIREBRICK 4');RGB = (/139,26,26/)
CASE ('FLESH');RGB = (/255,125,64/)
CASE ('FLORAL WHITE');RGB = (/255,250,240/)
CASE ('FOREST GREEN');RGB = (/34,139,34/)
CASE ('GAINSBORO');RGB = (/220,220,220/)
CASE ('GHOST WHITE');RGB = (/248,248,255/)
CASE ('GOLD');RGB = (/255,215,0/)
CASE ('GOLD 1');RGB = (/238,201,0/)
CASE ('GOLD 2');RGB = (/205,173,0/)
CASE ('GOLD 3');RGB = (/139,117,0/)
CASE ('GOLDENROD');RGB = (/218,165,32/)
CASE ('GOLDENROD 1');RGB = (/255,193,37/)
CASE ('GOLDENROD 2');RGB = (/238,180,34/)
CASE ('GOLDENROD 3');RGB = (/205,155,29/)
CASE ('GOLDENROD 4');RGB = (/139,105,20/)
CASE ('GRAY');RGB = (/128,128,128/)
CASE ('GRAY 1');RGB = (/3,3,3/)
CASE ('GRAY 10');RGB = (/26,26,26/)
CASE ('GRAY 11');RGB = (/28,28,28/)
CASE ('GRAY 12');RGB = (/31,31,31/)
CASE ('GRAY 13');RGB = (/33,33,33/)
CASE ('GRAY 14');RGB = (/36,36,36/)
CASE ('GRAY 15');RGB = (/38,38,38/)
CASE ('GRAY 16');RGB = (/41,41,41/)
CASE ('GRAY 17');RGB = (/43,43,43/)
CASE ('GRAY 18');RGB = (/46,46,46/)
CASE ('GRAY 19');RGB = (/48,48,48/)
CASE ('GRAY 2');RGB = (/5,5,5/)
CASE ('GRAY 20');RGB = (/51,51,51/)
CASE ('GRAY 21');RGB = (/54,54,54/)
CASE ('GRAY 22');RGB = (/56,56,56/)
CASE ('GRAY 23');RGB = (/59,59,59/)
CASE ('GRAY 24');RGB = (/61,61,61/)
CASE ('GRAY 25');RGB = (/64,64,64/)
CASE ('GRAY 26');RGB = (/66,66,66/)
CASE ('GRAY 27');RGB = (/69,69,69/)
CASE ('GRAY 28');RGB = (/71,71,71/)
CASE ('GRAY 29');RGB = (/74,74,74/)
CASE ('GRAY 3');RGB = (/8,8,8/)
CASE ('GRAY 30');RGB = (/77,77,77/)
CASE ('GRAY 31');RGB = (/79,79,79/)
CASE ('GRAY 32');RGB = (/82,82,82/)
CASE ('GRAY 33');RGB = (/84,84,84/)
CASE ('GRAY 34');RGB = (/87,87,87/)
CASE ('GRAY 35');RGB = (/89,89,89/)
CASE ('GRAY 36');RGB = (/92,92,92/)
CASE ('GRAY 37');RGB = (/94,94,94/)
CASE ('GRAY 38');RGB = (/97,97,97/)
CASE ('GRAY 39');RGB = (/99,99,99/)
CASE ('GRAY 4');RGB = (/10,10,10/)
CASE ('GRAY 40');RGB = (/102,102,102/)
CASE ('GRAY 42');RGB = (/107,107,107/)
CASE ('GRAY 43');RGB = (/110,110,110/)
CASE ('GRAY 44');RGB = (/112,112,112/)
CASE ('GRAY 45');RGB = (/115,115,115/)
CASE ('GRAY 46');RGB = (/117,117,117/)
CASE ('GRAY 47');RGB = (/120,120,120/)
CASE ('GRAY 48');RGB = (/122,122,122/)
CASE ('GRAY 49');RGB = (/125,125,125/)
CASE ('GRAY 5');RGB = (/13,13,13/)
CASE ('GRAY 50');RGB = (/127,127,127/)
CASE ('GRAY 51');RGB = (/130,130,130/)
CASE ('GRAY 52');RGB = (/133,133,133/)
CASE ('GRAY 53');RGB = (/135,135,135/)
CASE ('GRAY 54');RGB = (/138,138,138/)
CASE ('GRAY 55');RGB = (/140,140,140/)
CASE ('GRAY 56');RGB = (/143,143,143/)
CASE ('GRAY 57');RGB = (/145,145,145/)
CASE ('GRAY 58');RGB = (/148,148,148/)
CASE ('GRAY 59');RGB = (/150,150,150/)
CASE ('GRAY 6');RGB = (/15,15,15/)
CASE ('GRAY 60');RGB = (/153,153,153/)
CASE ('GRAY 61');RGB = (/156,156,156/)
CASE ('GRAY 62');RGB = (/158,158,158/)
CASE ('GRAY 63');RGB = (/161,161,161/)
CASE ('GRAY 64');RGB = (/163,163,163/)
CASE ('GRAY 65');RGB = (/166,166,166/)
CASE ('GRAY 66');RGB = (/168,168,168/)
CASE ('GRAY 67');RGB = (/171,171,171/)
CASE ('GRAY 68');RGB = (/173,173,173/)
CASE ('GRAY 69');RGB = (/176,176,176/)
CASE ('GRAY 7');RGB = (/18,18,18/)
CASE ('GRAY 70');RGB = (/179,179,179/)
CASE ('GRAY 71');RGB = (/181,181,181/)
CASE ('GRAY 72');RGB = (/184,184,184/)
CASE ('GRAY 73');RGB = (/186,186,186/)
CASE ('GRAY 74');RGB = (/189,189,189/)
CASE ('GRAY 75');RGB = (/191,191,191/)
CASE ('GRAY 76');RGB = (/194,194,194/)
CASE ('GRAY 77');RGB = (/196,196,196/)
CASE ('GRAY 78');RGB = (/199,199,199/)
CASE ('GRAY 79');RGB = (/201,201,201/)
CASE ('GRAY 8');RGB = (/20,20,20/)
CASE ('GRAY 80');RGB = (/204,204,204/)
CASE ('GRAY 81');RGB = (/207,207,207/)
CASE ('GRAY 82');RGB = (/209,209,209/)
CASE ('GRAY 83');RGB = (/212,212,212/)
CASE ('GRAY 84');RGB = (/214,214,214/)
CASE ('GRAY 85');RGB = (/217,217,217/)
CASE ('GRAY 86');RGB = (/219,219,219/)
CASE ('GRAY 87');RGB = (/222,222,222/)
CASE ('GRAY 88');RGB = (/224,224,224/)
CASE ('GRAY 89');RGB = (/227,227,227/)
CASE ('GRAY 9');RGB = (/23,23,23/)
CASE ('GRAY 90');RGB = (/229,229,229/)
CASE ('GRAY 91');RGB = (/232,232,232/)
CASE ('GRAY 92');RGB = (/235,235,235/)
CASE ('GRAY 93');RGB = (/237,237,237/)
CASE ('GRAY 94');RGB = (/240,240,240/)
CASE ('GRAY 95');RGB = (/242,242,242/)
CASE ('GRAY 97');RGB = (/247,247,247/)
CASE ('GRAY 98');RGB = (/250,250,250/)
CASE ('GRAY 99');RGB = (/252,252,252/)
CASE ('GREEN');RGB = (/0,255,0/)
CASE ('GREEN 2');RGB = (/0,238,0/)
CASE ('GREEN 3');RGB = (/0,205,0/)
CASE ('GREEN 4');RGB = (/0,139,0/)
CASE ('GREEN YELLOW');RGB = (/173,255,47/)
CASE ('HONEYDEW');RGB = (/240,255,240/)
CASE ('HONEYDEW 1');RGB = (/224,238,224/)
CASE ('HONEYDEW 2');RGB = (/193,205,193/)
CASE ('HONEYDEW 3');RGB = (/131,139,131/)
CASE ('HOT PINK');RGB = (/255,105,180/)
CASE ('HOT PINK 1');RGB = (/255,110,180/)
CASE ('HOT PINK 2');RGB = (/238,106,167/)
CASE ('HOT PINK 3');RGB = (/205,96,144/)
CASE ('HOT PINK 4');RGB = (/139,58,98/)
CASE ('INDIAN RED');RGB = (/205,92,92/)
CASE ('INDIAN RED 1');RGB = (/255,106,106/)
CASE ('INDIAN RED 2');RGB = (/238,99,99/)
CASE ('INDIAN RED 3');RGB = (/205,85,85/)
CASE ('INDIAN RED 4');RGB = (/139,58,58/)
CASE ('INDIGO');RGB = (/75,0,130/)
CASE ('IVORY');RGB = (/255,255,240/)
CASE ('IVORY 1');RGB = (/238,238,224/)
CASE ('IVORY 2');RGB = (/205,205,193/)
CASE ('IVORY 3');RGB = (/139,139,131/)
CASE ('IVORY BLACK');RGB = (/41,36,33/)
CASE ('KELLY GREEN');RGB = (/0,128,0/)
CASE ('KHAKI');RGB = (/240,230,140/)
CASE ('KHAKI 1');RGB = (/255,246,143/)
CASE ('KHAKI 2');RGB = (/238,230,133/)
CASE ('KHAKI 3');RGB = (/205,198,115/)
CASE ('KHAKI 4');RGB = (/139,134,78/)
CASE ('LAVENDER');RGB = (/230,230,250/)
CASE ('LAVENDER BLUSH');RGB = (/255,240,245/)
CASE ('LAVENDER BLUSH 1');RGB = (/238,224,229/)
CASE ('LAVENDER BLUSH 2');RGB = (/205,193,197/)
CASE ('LAVENDER BLUSH 3');RGB = (/139,131,134/)
CASE ('LAWN GREEN');RGB = (/124,252,0/)
CASE ('LEMON CHIFFON');RGB = (/255,250,205/)
CASE ('LEMON CHIFFON 1');RGB = (/238,233,191/)
CASE ('LEMON CHIFFON 2');RGB = (/205,201,165/)
CASE ('LEMON CHIFFON 3');RGB = (/139,137,112/)
CASE ('LIGHT BLUE');RGB = (/173,216,230/)
CASE ('LIGHT BLUE 1');RGB = (/191,239,255/)
CASE ('LIGHT BLUE 2');RGB = (/178,223,238/)
CASE ('LIGHT BLUE 3');RGB = (/154,192,205/)
CASE ('LIGHT BLUE 4');RGB = (/104,131,139/)
CASE ('LIGHT CORAL');RGB = (/240,128,128/)
CASE ('LIGHT CYAN');RGB = (/224,255,255/)
CASE ('LIGHT CYAN 1');RGB = (/209,238,238/)
CASE ('LIGHT CYAN 2');RGB = (/180,205,205/)
CASE ('LIGHT CYAN 3');RGB = (/122,139,139/)
CASE ('LIGHT GOLDENROD');RGB = (/255,236,139/)
CASE ('LIGHT GOLDENROD 1');RGB = (/238,220,130/)
CASE ('LIGHT GOLDENROD 2');RGB = (/205,190,112/)
CASE ('LIGHT GOLDENROD 3');RGB = (/139,129,76/)
CASE ('LIGHT GOLDENROD YELLOW');RGB = (/250,250,210/)
CASE ('LIGHT GREY');RGB = (/211,211,211/)
CASE ('LIGHT PINK');RGB = (/255,182,193/)
CASE ('LIGHT PINK 1');RGB = (/255,174,185/)
CASE ('LIGHT PINK 2');RGB = (/238,162,173/)
CASE ('LIGHT PINK 3');RGB = (/205,140,149/)
CASE ('LIGHT PINK 4');RGB = (/139,95,101/)
CASE ('LIGHT SALMON');RGB = (/255,160,122/)
CASE ('LIGHT SALMON 1');RGB = (/238,149,114/)
CASE ('LIGHT SALMON 2');RGB = (/205,129,98/)
CASE ('LIGHT SALMON 3');RGB = (/139,87,66/)
CASE ('LIGHT SEA GREEN');RGB = (/32,178,170/)
CASE ('LIGHT SKY BLUE');RGB = (/135,206,250/)
CASE ('LIGHT SKY BLUE 1');RGB = (/176,226,255/)
CASE ('LIGHT SKY BLUE 2');RGB = (/164,211,238/)
CASE ('LIGHT SKY BLUE 3');RGB = (/141,182,205/)
CASE ('LIGHT SKY BLUE 4');RGB = (/96,123,139/)
CASE ('LIGHT SLATE BLUE');RGB = (/132,112,255/)
CASE ('LIGHT SLATE GRAY');RGB = (/119,136,153/)
CASE ('LIGHT STEEL BLUE');RGB = (/176,196,222/)
CASE ('LIGHT STEEL BLUE 1');RGB = (/202,225,255/)
CASE ('LIGHT STEEL BLUE 2');RGB = (/188,210,238/)
CASE ('LIGHT STEEL BLUE 3');RGB = (/162,181,205/)
CASE ('LIGHT STEEL BLUE 4');RGB = (/110,123,139/)
CASE ('LIGHT YELLOW 1');RGB = (/255,255,224/)
CASE ('LIGHT YELLOW 2');RGB = (/238,238,209/)
CASE ('LIGHT YELLOW 3');RGB = (/205,205,180/)
CASE ('LIGHT YELLOW 4');RGB = (/139,139,122/)
CASE ('LIME GREEN');RGB = (/50,205,50/)
CASE ('LINEN');RGB = (/250,240,230/)
CASE ('MAGENTA');RGB = (/255,0,255/)
CASE ('MAGENTA 2');RGB = (/238,0,238/)
CASE ('MAGENTA 3');RGB = (/205,0,205/)
CASE ('MAGENTA 4');RGB = (/139,0,139/)
CASE ('MANGANESE BLUE');RGB = (/3,168,158/)
CASE ('MAROON');RGB = (/128,0,0/)
CASE ('MAROON 1');RGB = (/255,52,179/)
CASE ('MAROON 2');RGB = (/238,48,167/)
CASE ('MAROON 3');RGB = (/205,41,144/)
CASE ('MAROON 4');RGB = (/139,28,98/)
CASE ('MEDIUM ORCHID');RGB = (/186,85,211/)
CASE ('MEDIUM ORCHID 1');RGB = (/224,102,255/)
CASE ('MEDIUM ORCHID 2');RGB = (/209,95,238/)
CASE ('MEDIUM ORCHID 3');RGB = (/180,82,205/)
CASE ('MEDIUM ORCHID 4');RGB = (/122,55,139/)
CASE ('MEDIUM PURPLE');RGB = (/147,112,219/)
CASE ('MEDIUM PURPLE 1');RGB = (/171,130,255/)
CASE ('MEDIUM PURPLE 2');RGB = (/159,121,238/)
CASE ('MEDIUM PURPLE 3');RGB = (/137,104,205/)
CASE ('MEDIUM PURPLE 4');RGB = (/93,71,139/)
CASE ('MEDIUM SEA GREEN');RGB = (/60,179,113/)
CASE ('MEDIUM SLATE BLUE');RGB = (/123,104,238/)
CASE ('MEDIUM SPRING GREEN');RGB = (/0,250,154/)
CASE ('MEDIUM TURQUOISE');RGB = (/72,209,204/)
CASE ('MEDIUM VIOLET RED');RGB = (/199,21,133/)
CASE ('MELON');RGB = (/227,168,105/)
CASE ('MIDNIGHT BLUE');RGB = (/25,25,112/)
CASE ('MINT');RGB = (/189,252,201/)
CASE ('MINT CREAM');RGB = (/245,255,250/)
CASE ('MISTY ROSE');RGB = (/255,228,225/)
CASE ('MISTY ROSE 1');RGB = (/238,213,210/)
CASE ('MISTY ROSE 2');RGB = (/205,183,181/)
CASE ('MISTY ROSE 3');RGB = (/139,125,123/)
CASE ('MOCCASIN');RGB = (/255,228,181/)
CASE ('NAVAJO WHITE');RGB = (/255,222,173/)
CASE ('NAVAJO WHITE 1');RGB = (/238,207,161/)
CASE ('NAVAJO WHITE 2');RGB = (/205,179,139/)
CASE ('NAVAJO WHITE 3');RGB = (/139,121,94/)
CASE ('NAVY');RGB = (/0,0,128/)
CASE ('OLD LACE');RGB = (/253,245,230/)
CASE ('OLIVE');RGB = (/128,128,0/)
CASE ('OLIVE DRAB');RGB = (/192,255,62/)
CASE ('OLIVE DRAB 1');RGB = (/179,238,58/)
CASE ('OLIVE DRAB 2');RGB = (/154,205,50/)
CASE ('OLIVE DRAB 3');RGB = (/105,139,34/)
CASE ('ORANGE');RGB = (/255,128,0/)
CASE ('ORANGE 1');RGB = (/255,165,0/)
CASE ('ORANGE 2');RGB = (/238,154,0/)
CASE ('ORANGE 3');RGB = (/205,133,0/)
CASE ('ORANGE 4');RGB = (/139,90,0/)
CASE ('ORANGE RED');RGB = (/255,69,0/)
CASE ('ORANGE RED 1');RGB = (/238,64,0/)
CASE ('ORANGE RED 2');RGB = (/205,55,0/)
CASE ('ORANGE RED 3');RGB = (/139,37,0/)
CASE ('ORCHID');RGB = (/218,112,214/)
CASE ('ORCHID 1');RGB = (/255,131,250/)
CASE ('ORCHID 2');RGB = (/238,122,233/)
CASE ('ORCHID 3');RGB = (/205,105,201/)
CASE ('ORCHID 4');RGB = (/139,71,137/)
CASE ('PALE GOLDENROD');RGB = (/238,232,170/)
CASE ('PALE GREEN');RGB = (/152,251,152/)
CASE ('PALE GREEN 1');RGB = (/154,255,154/)
CASE ('PALE GREEN 2');RGB = (/144,238,144/)
CASE ('PALE GREEN 3');RGB = (/124,205,124/)
CASE ('PALE GREEN 4');RGB = (/84,139,84/)
CASE ('PALE TURQUOISE');RGB = (/187,255,255/)
CASE ('PALE TURQUOISE 1');RGB = (/174,238,238/)
CASE ('PALE TURQUOISE 2');RGB = (/150,205,205/)
CASE ('PALE TURQUOISE 3');RGB = (/102,139,139/)
CASE ('PALE VIOLET RED');RGB = (/219,112,147/)
CASE ('PALE VIOLET RED 1');RGB = (/255,130,171/)
CASE ('PALE VIOLET RED 2');RGB = (/238,121,159/)
CASE ('PALE VIOLET RED 3');RGB = (/205,104,137/)
CASE ('PALE VIOLET RED 4');RGB = (/139,71,93/)
CASE ('PAPAYA WHIP');RGB = (/255,239,213/)
CASE ('PEACH PUFF');RGB = (/255,218,185/)
CASE ('PEACH PUFF 1');RGB = (/238,203,173/)
CASE ('PEACH PUFF 2');RGB = (/205,175,149/)
CASE ('PEACH PUFF 3');RGB = (/139,119,101/)
CASE ('PEACOCK');RGB = (/51,161,201/)
CASE ('PINK');RGB = (/255,192,203/)
CASE ('PINK 1');RGB = (/255,181,197/)
CASE ('PINK 2');RGB = (/238,169,184/)
CASE ('PINK 3');RGB = (/205,145,158/)
CASE ('PINK 4');RGB = (/139,99,108/)
CASE ('PLUM');RGB = (/221,160,221/)
CASE ('PLUM 1');RGB = (/255,187,255/)
CASE ('PLUM 2');RGB = (/238,174,238/)
CASE ('PLUM 3');RGB = (/205,150,205/)
CASE ('PLUM 4');RGB = (/139,102,139/)
CASE ('POWDER BLUE');RGB = (/176,224,230/)
CASE ('PURPLE');RGB = (/128,0,128/)
CASE ('PURPLE 1');RGB = (/155,48,255/)
CASE ('PURPLE 2');RGB = (/145,44,238/)
CASE ('PURPLE 3');RGB = (/125,38,205/)
CASE ('PURPLE 4');RGB = (/85,26,139/)
CASE ('RASPBERRY');RGB = (/135,38,87/)
CASE ('RAW SIENNA');RGB = (/199,97,20/)
CASE ('RED');RGB = (/255,0,0/)
CASE ('RED 1');RGB = (/238,0,0/)
CASE ('RED 2');RGB = (/205,0,0/)
CASE ('RED 3');RGB = (/139,0,0/)
CASE ('ROSY BROWN');RGB = (/188,143,143/)
CASE ('ROSY BROWN 1');RGB = (/255,193,193/)
CASE ('ROSY BROWN 2');RGB = (/238,180,180/)
CASE ('ROSY BROWN 3');RGB = (/205,155,155/)
CASE ('ROSY BROWN 4');RGB = (/139,105,105/)
CASE ('ROYAL BLUE');RGB = (/65,105,225/)
CASE ('ROYAL BLUE 1');RGB = (/72,118,255/)
CASE ('ROYAL BLUE 2');RGB = (/67,110,238/)
CASE ('ROYAL BLUE 3');RGB = (/58,95,205/)
CASE ('ROYAL BLUE 4');RGB = (/39,64,139/)
CASE ('SALMON');RGB = (/250,128,114/)
CASE ('SALMON 1');RGB = (/255,140,105/)
CASE ('SALMON 2');RGB = (/238,130,98/)
CASE ('SALMON 3');RGB = (/205,112,84/)
CASE ('SALMON 4');RGB = (/139,76,57/)
CASE ('SANDY BROWN');RGB = (/244,164,96/)
CASE ('SAP GREEN');RGB = (/48,128,20/)
CASE ('SEA GREEN');RGB = (/84,255,159/)
CASE ('SEA GREEN 1');RGB = (/78,238,148/)
CASE ('SEA GREEN 2');RGB = (/67,205,128/)
CASE ('SEA GREEN 3');RGB = (/46,139,87/)
CASE ('SEASHELL');RGB = (/255,245,238/)
CASE ('SEASHELL 1');RGB = (/238,229,222/)
CASE ('SEASHELL 2');RGB = (/205,197,191/)
CASE ('SEASHELL 3');RGB = (/139,134,130/)
CASE ('SEPIA');RGB = (/94,38,18/)
CASE ('SIENNA');RGB = (/160,82,45/)
CASE ('SIENNA 1');RGB = (/255,130,71/)
CASE ('SIENNA 2');RGB = (/238,121,66/)
CASE ('SIENNA 3');RGB = (/205,104,57/)
CASE ('SIENNA 4');RGB = (/139,71,38/)
CASE ('SILVER');RGB = (/192,192,192/)
CASE ('SKY BLUE');RGB = (/135,206,235/)
CASE ('SKY BLUE 1');RGB = (/135,206,255/)
CASE ('SKY BLUE 2');RGB = (/126,192,238/)
CASE ('SKY BLUE 3');RGB = (/108,166,205/)
CASE ('SKY BLUE 4');RGB = (/74,112,139/)
CASE ('SLATE BLUE');RGB = (/106,90,205/)
CASE ('SLATE BLUE 1');RGB = (/131,111,255/)
CASE ('SLATE BLUE 2');RGB = (/122,103,238/)
CASE ('SLATE BLUE 3');RGB = (/105,89,205/)
CASE ('SLATE BLUE 4');RGB = (/71,60,139/)
CASE ('SLATE GRAY');RGB = (/112,128,144/)
CASE ('SLATE GRAY 1');RGB = (/198,226,255/)
CASE ('SLATE GRAY 2');RGB = (/185,211,238/)
CASE ('SLATE GRAY 3');RGB = (/159,182,205/)
CASE ('SLATE GRAY 4');RGB = (/108,123,139/)
CASE ('SNOW');RGB = (/255,250,250/)
CASE ('SNOW 1');RGB = (/238,233,233/)
CASE ('SNOW 2');RGB = (/205,201,201/)
CASE ('SNOW 3');RGB = (/139,137,137/)
CASE ('SPRING GREEN');RGB = (/0,255,127/)
CASE ('SPRING GREEN 1');RGB = (/0,238,118/)
CASE ('SPRING GREEN 2');RGB = (/0,205,102/)
CASE ('SPRING GREEN 3');RGB = (/0,139,69/)
CASE ('STEEL BLUE');RGB = (/70,130,180/)
CASE ('STEEL BLUE 1');RGB = (/99,184,255/)
CASE ('STEEL BLUE 2');RGB = (/92,172,238/)
CASE ('STEEL BLUE 3');RGB = (/79,148,205/)
CASE ('STEEL BLUE 4');RGB = (/54,100,139/)
CASE ('TAN');RGB = (/210,180,140/)
CASE ('TAN 1');RGB = (/255,165,79/)
CASE ('TAN 2');RGB = (/238,154,73/)
CASE ('TAN 3');RGB = (/205,133,63/)
CASE ('TAN 4');RGB = (/139,90,43/)
CASE ('TEAL');RGB = (/0,128,128/)
CASE ('THISTLE');RGB = (/216,191,216/)
CASE ('THISTLE 1');RGB = (/255,225,255/)
CASE ('THISTLE 2');RGB = (/238,210,238/)
CASE ('THISTLE 3');RGB = (/205,181,205/)
CASE ('THISTLE 4');RGB = (/139,123,139/)
CASE ('TOMATO');RGB = (/255,99,71/)
CASE ('TOMATO 1');RGB = (/238,92,66/)
CASE ('TOMATO 2');RGB = (/205,79,57/)
CASE ('TOMATO 3');RGB = (/139,54,38/)
CASE ('TURQUOISE');RGB = (/64,224,208/)
CASE ('TURQUOISE 1');RGB = (/0,245,255/)
CASE ('TURQUOISE 2');RGB = (/0,229,238/)
CASE ('TURQUOISE 3');RGB = (/0,197,205/)
CASE ('TURQUOISE 4');RGB = (/0,134,139/)
CASE ('TURQUOISE BLUE');RGB = (/0,199,140/)
CASE ('VIOLET');RGB = (/238,130,238/)
CASE ('VIOLET RED');RGB = (/208,32,144/)
CASE ('VIOLET RED 1');RGB = (/255,62,150/)
CASE ('VIOLET RED 2');RGB = (/238,58,140/)
CASE ('VIOLET RED 3');RGB = (/205,50,120/)
CASE ('VIOLET RED 4');RGB = (/139,34,82/)
CASE ('WARM GREY');RGB = (/128,128,105/)
CASE ('WHEAT');RGB = (/245,222,179/)
CASE ('WHEAT 1');RGB = (/255,231,186/)
CASE ('WHEAT 2');RGB = (/238,216,174/)
CASE ('WHEAT 3');RGB = (/205,186,150/)
CASE ('WHEAT 4');RGB = (/139,126,102/)
CASE ('WHITE');RGB = (/255,255,255/)
CASE ('WHITE SMOKE');RGB = (/245,245,245/)
CASE ('YELLOW');RGB = (/255,255,0/)
CASE ('YELLOW 1');RGB = (/238,238,0/)
CASE ('YELLOW 2');RGB = (/205,205,0/)
CASE ('YELLOW 3');RGB = (/139,139,0/)

CASE DEFAULT
   WRITE(MESSAGE,'(A,A,A)') "ERROR: The COLOR, ", TRIM(COLOR),", is not a defined color"
   CALL SHUTDOWN(MESSAGE)      
END SELECT

END SUBROUTINE COLOR2RGB


SUBROUTINE CHANGE_UNITS(QUANTITY,UNITS,STATISTICS,MYID,LU_ERR)
CHARACTER(30), INTENT(IN) :: QUANTITY,STATISTICS
INTEGER, INTENT(IN) :: MYID,LU_ERR
CHARACTER(30), INTENT(INOUT) ::UNITS
CHARACTER(150) :: MESSAGE
INTEGER :: I,UNIT_L

UNIT_L = LEN(TRIM(UNITS))

SELECT CASE (STATISTICS)
   CASE('VOLUME INTEGRAL')
      I = INDEX(UNITS,'/m3')
      IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
   CASE('AREA INTEGRAL')
      I = INDEX(UNITS,'/m2')
      IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
   CASE('SURFACE INTEGRAL')
      I = INDEX(UNITS,'/m2')
      IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
   CASE('TIME INTEGRAL')
      I = INDEX(UNITS,'/s')
      IF (I/=0) THEN
         WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+2:UNIT_L)
      ELSE
         I = INDEX(UNITS,'kW')
         IF (I==1) WRITE(UNITS,'(A,A)') 'kJ',UNITS(3:UNIT_L)
         IF (I>1)  WRITE(UNITS,'(A,A,A)') UNITS(1:I-1),'kJ',UNITS(I+2:UNIT_L)
      ENDIF
      IF (QUANTITY=='FED') I=1
   CASE DEFAULT
      RETURN
END SELECT

IF (I==0) THEN
   WRITE(MESSAGE,'(A,A,A,A)')  'WARNING: Problem with units compatibility using STATISTICS of ',TRIM(STATISTICS), &
                               ' with the QUANTITY ',TRIM(QUANTITY)
   IF (MYID==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE)
ENDIF

RETURN

END SUBROUTINE CHANGE_UNITS


END MODULE COMP_FUNCTIONS



MODULE MEMORY_FUNCTIONS

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE COMP_FUNCTIONS, ONLY : SHUTDOWN
IMPLICIT NONE

CONTAINS


SUBROUTINE ChkMemErr(CodeSect,VarName,IZERO)
 
! Memory checking routine
 
CHARACTER(*), INTENT(IN) :: CodeSect, VarName
INTEGER IZERO
CHARACTER(100) MESSAGE
 
IF (IZERO==0) RETURN
 
WRITE(MESSAGE,'(4A)') 'ERROR: Memory allocation failed for ', TRIM(VarName),' in the routine ',TRIM(CodeSect)
CALL SHUTDOWN(MESSAGE)

END SUBROUTINE ChkMemErr



SUBROUTINE RE_ALLOCATE_DROPLETS(CODE,NM,NOM,NEW_DROPS)
 
TYPE (DROPLET_TYPE), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: CODE,NM,NOM,NEW_DROPS
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE(OMESH_TYPE), POINTER :: M2=>NULL()
 
SELECT CASE(CODE)
   CASE(1)
      M=>MESHES(NM)
      ALLOCATE(DUMMY(1:M%NLPDIM))
      DUMMY = M%DROPLET
      DEALLOCATE(M%DROPLET)
      ALLOCATE(M%DROPLET(M%NLPDIM+NEW_DROPS))
      M%DROPLET(1:M%NLPDIM) = DUMMY(1:M%NLPDIM)
      M%NLPDIM = M%NLPDIM+NEW_DROPS
   CASE(2)
      M2=>MESHES(NM)%OMESH(NOM)
      ALLOCATE(DUMMY(1:M2%N_DROP_ORPHANS_DIM))
      DUMMY = M2%DROPLET
      DEALLOCATE(M2%DROPLET)
      ALLOCATE(M2%DROPLET(M2%N_DROP_ORPHANS_DIM+NEW_DROPS))
      M2%DROPLET(1:M2%N_DROP_ORPHANS_DIM) = DUMMY(1:M2%N_DROP_ORPHANS_DIM)
      M2%N_DROP_ORPHANS_DIM = M2%N_DROP_ORPHANS_DIM + NEW_DROPS
END SELECT
DEALLOCATE(DUMMY)

END SUBROUTINE RE_ALLOCATE_DROPLETS

 

FUNCTION REALLOCATE(P,N1,N2)          
REAL(EB), POINTER, DIMENSION(:) :: P, REALLOCATE
INTEGER, INTENT(IN) :: N1,N2
INTEGER :: NOLD, IERR
CHARACTER(100) :: MESSAGE
ALLOCATE(REALLOCATE(N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE'
   CALL SHUTDOWN(MESSAGE)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
NOLD = MIN(SIZE(P), N2-N1+1)
REALLOCATE(N1:NOLD+N1-1) = P(N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P) 
END FUNCTION REALLOCATE



SUBROUTINE RE_ALLOCATE_STRINGS(NM)
 
CHARACTER(80), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M=>NULL()
 
M=>MESHES(NM)
ALLOCATE(DUMMY(1:M%N_STRINGS))
DUMMY = M%STRING
DEALLOCATE(M%STRING)
ALLOCATE(M%STRING(M%N_STRINGS_MAX+100))
M%STRING(1:M%N_STRINGS) = DUMMY(1:M%N_STRINGS)
M%N_STRINGS_MAX = M%N_STRINGS_MAX+100
DEALLOCATE(DUMMY)
 
END SUBROUTINE RE_ALLOCATE_STRINGS



SUBROUTINE ALLOCATE_WALL_ARRAYS(NM,IW)

USE GLOBAL_CONSTANTS, ONLY: RADIATION,SIGMA,TMPA4,OPEN_BOUNDARY,EVACUATION_ONLY,&
                            NRA=>NUMBER_RADIATION_ANGLES,NSB=>NUMBER_SPECTRAL_BANDS, &
                            TMPA
INTEGER, INTENT(IN) :: NM,IW
INTEGER :: IZERO,II,IL,NN,N
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE (SURFACE_TYPE), POINTER :: SF=>NULL()
TYPE (WALL_TYPE), POINTER :: WC=>NULL()

M=>MESHES(NM)
SF => SURFACE(M%IJKW(5,IW))
WC => M%WALL(IW)
IF (EVACUATION_ONLY(NM)) RETURN
IF (M%BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) RETURN

IF (RADIATION) THEN
   ALLOCATE(WC%ILW(NRA,NSB),STAT=IZERO)
   CALL ChkMemErr('FUNC','ILW',IZERO)
   WC%ILW        = SIGMA*TMPA4*RPI
   M%QRADIN(IW)  = M%E_WALL(IW)*SIGMA*TMPA4
   M%QRADOUT(IW) = M%E_WALL(IW)*SIGMA*TMPA4
ENDIF

IF (SF%THERMALLY_THICK) THEN
   ALLOCATE(WC%TMP_S(0:SF%N_CELLS+1),STAT=IZERO)
   CALL ChkMemErr('FUNC','TMP_S',IZERO)
   WC%TMP_S = SF%TMP_INNER
   ALLOCATE(WC%RHO_S(0:SF%N_CELLS+1,SF%N_MATL),STAT=IZERO)
   CALL ChkMemErr('FUNC','RHO_S',IZERO)
   WC%RHO_S = 0._EB
   DO II=0,SF%N_CELLS+1
      IL = SF%LAYER_INDEX(II)
      DO NN=1,SF%N_LAYER_MATL(IL)
         DO N=1,SF%N_MATL
            IF (SF%LAYER_MATL_INDEX(IL,NN)==SF%MATL_INDEX(N)) &
               WC%RHO_S(II,N) = SF%LAYER_MATL_FRAC(IL,NN)*SF%LAYER_DENSITY(IL)
         ENDDO
      ENDDO
   ENDDO
   IF (SF%SHRINK) THEN
      ALLOCATE(WC%LAYER_THICKNESS(SF%N_LAYERS),STAT=IZERO)
      CALL ChkMemErr('FUNC','LAYER_THICKNESS',IZERO)
      WC%LAYER_THICKNESS = SF%LAYER_THICKNESS(1:SF%N_LAYERS)
      ALLOCATE(WC%N_LAYER_CELLS(SF%N_LAYERS),STAT=IZERO)
      CALL ChkMemErr('FUNC','N_LAYER_CELLS',IZERO)
      WC%N_LAYER_CELLS = SF%N_LAYER_CELLS
      ALLOCATE(WC%X_S(0:SF%N_CELLS),STAT=IZERO)
      CALL ChkMemErr('FUNC','X_S',IZERO)
      WC%X_S = SF%X_S
   ENDIF
   WC%SHRINKING = .FALSE.
   WC%BURNAWAY  = .FALSE.
   M%TMP_F(IW) = SF%TMP_INNER
   M%TMP_B(IW) = SF%TMP_INNER
ENDIF

IF (SF%VEGETATION) THEN

   ALLOCATE(WC%VEG_FUELMASS_L(SF%NVEG_L),STAT=IZERO)
   CALL ChkMemErr('FUNC','VEG_FUELMASS_L',IZERO)
   ALLOCATE(WC%VEG_MOISTMASS_L(SF%NVEG_L),STAT=IZERO)
   CALL ChkMemErr('FUNC','VEG_MOISTMASS_L',IZERO)
   ALLOCATE(WC%VEG_TMP_L(0:SF%NVEG_L),STAT=IZERO)
   CALL ChkMemErr('FUNC','VEG_TMP_L',IZERO)

   WC%VEG_FUELMASS_L(:)  = SF%VEG_LOAD / SF%NVEG_L
   WC%VEG_MOISTMASS_L(:) = SF%VEG_MOISTURE*WC%VEG_FUELMASS_L(:)
   WC%VEG_TMP_L(:)       = SF%VEG_INITIAL_TEMP + 273.15_EB
   M%TMP_F(IW)           = SF%VEG_INITIAL_TEMP + 273.15_EB
   IF(SF%VEG_NO_BURN) M%TMP_F(IW) = TMPA

ENDIF

WC%ALREADY_ALLOCATED = .TRUE.

END SUBROUTINE ALLOCATE_WALL_ARRAYS

END MODULE MEMORY_FUNCTIONS 



MODULE GEOMETRY_FUNCTIONS

! Functions for manipulating geometry

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
IMPLICIT NONE

CONTAINS
 

SUBROUTINE SEARCH_OTHER_MESHES(XX,YY,ZZ,NOM,IIO,JJO,KKO)

! Given the point (XX,YY,ZZ), determine which other mesh it intersects and what its indices are.

REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(OUT) :: NOM,IIO,JJO,KKO
TYPE (MESH_TYPE), POINTER :: M2=>NULL()

OTHER_MESH_LOOP: DO NOM=1,NMESHES
   IF (EVACUATION_ONLY(NOM)) CYCLE OTHER_MESH_LOOP
   M2=>MESHES(NOM)
   IF (XX>=M2%XS .AND. XX<=M2%XF .AND.  YY>=M2%YS .AND. YY<=M2%YF .AND. ZZ>=M2%ZS .AND. ZZ<=M2%ZF) THEN
      XI  = MAX( 1._EB , MIN( REAL(M2%IBAR,EB)+ALMOST_ONE , M2%CELLSI(FLOOR((XX-M2%XS)*M2%RDXINT)) + 1._EB ) )
      YJ  = MAX( 1._EB , MIN( REAL(M2%JBAR,EB)+ALMOST_ONE , M2%CELLSJ(FLOOR((YY-M2%YS)*M2%RDYINT)) + 1._EB ) )
      ZK  = MAX( 1._EB , MIN( REAL(M2%KBAR,EB)+ALMOST_ONE , M2%CELLSK(FLOOR((ZZ-M2%ZS)*M2%RDZINT)) + 1._EB ) )
      IIO = FLOOR(XI)
      JJO = FLOOR(YJ)
      KKO = FLOOR(ZK)
      RETURN
   ENDIF
ENDDO OTHER_MESH_LOOP

NOM = 0

END SUBROUTINE SEARCH_OTHER_MESHES



SUBROUTINE BLOCK_CELL(NM,I1,I2,J1,J2,K1,K2,IVAL,OBST_INDEX)

! Indicate which cells are blocked off
 
INTEGER :: NM,I1,I2,J1,J2,K1,K2,IVAL,I,J,K,OBST_INDEX,IC
TYPE (MESH_TYPE), POINTER :: M=>NULL()
 
M => MESHES(NM)
DO K=K1,K2
   DO J=J1,J2
      DO I=I1,I2
         IC = M%CELL_INDEX(I,J,K)
         SELECT CASE(IVAL)
            CASE(0) 
               M%SOLID(IC)        = .FALSE.
               M%OBST_INDEX_C(IC) = 0
            CASE(1)
               M%SOLID(IC)        = .TRUE.
               M%OBST_INDEX_C(IC) = OBST_INDEX
         END SELECT
         IF (OBST_INDEX==0) M%EXTERIOR(IC) = .TRUE.
      ENDDO
   ENDDO
ENDDO
 
END SUBROUTINE BLOCK_CELL
 

 
LOGICAL FUNCTION INTERIOR(XX,YY,ZZ)

INTEGER NM
REAL(EB), INTENT(IN) :: XX,YY,ZZ

INTERIOR = .FALSE.

DO NM=1,NMESHES
   IF (XX>MESHES(NM)%XS .AND. XX<MESHES(NM)%XF .AND. &
       YY>MESHES(NM)%YS .AND. YY<MESHES(NM)%YF .AND. &
       ZZ>MESHES(NM)%ZS .AND. ZZ<MESHES(NM)%ZF) INTERIOR = .TRUE.
ENDDO

END FUNCTION INTERIOR



SUBROUTINE ASSIGN_PRESSURE_ZONE(NM,XX,YY,ZZ,I_ZONE)

! Given the point (XX,YY,ZZ) within Mesh NM, determine all mesh cells within the same pressure zone

REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(IN) :: NM,I_ZONE
INTEGER :: NN,IOR,IC,II,JJ,KK,III,JJJ,KKK,Q_N,IIO,JJO,KKO,NOM
INTEGER, ALLOCATABLE, DIMENSION(:) :: Q_I,Q_J,Q_K
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE (OBSTRUCTION_TYPE), POINTER :: OB=>NULL()

IF (EVACUATION_ONLY(NM)) RETURN

M=>MESHES(NM)

ALLOCATE(Q_I(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_J(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_K(M%IBAR*M%JBAR*M%KBAR))

! Find the cell indices corresponding to the given point

XI  = MAX( 1._EB , MIN( REAL(M%IBAR,EB)+ALMOST_ONE , M%CELLSI(NINT((XX-M%XS)*M%RDXINT)) + 1._EB ) )
YJ  = MAX( 1._EB , MIN( REAL(M%JBAR,EB)+ALMOST_ONE , M%CELLSJ(NINT((YY-M%YS)*M%RDYINT)) + 1._EB ) )
ZK  = MAX( 1._EB , MIN( REAL(M%KBAR,EB)+ALMOST_ONE , M%CELLSK(NINT((ZZ-M%ZS)*M%RDZINT)) + 1._EB ) )
II  = FLOOR(XI)
JJ  = FLOOR(YJ)
KK  = FLOOR(ZK)

! Add the first entry to "queue" of cells that need a pressure zone number

Q_I(1) = II
Q_J(1) = JJ
Q_K(1) = KK
Q_N    = 1
M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE

! Look to all cells adjacent to the starting cell and determine if they are in the ZONE as well. 
! Repeat process until all cells are found.

SORT_QUEUE: DO 

   IF (Q_N<1) EXIT SORT_QUEUE

   III = Q_I(Q_N)
   JJJ = Q_J(Q_N)
   KKK = Q_K(Q_N)
   IC  = M%CELL_INDEX(III,JJJ,KKK)
   Q_N = Q_N - 1

   SEARCH_LOOP: DO IOR=-3,3

      IF (IOR==0) CYCLE SEARCH_LOOP

      SELECT CASE(IOR)
         CASE(-1)
            II = III-1
            JJ = JJJ
            KK = KKK
         CASE( 1)
            II = III+1
            JJ = JJJ
            KK = KKK
         CASE(-2)
            II = III
            JJ = JJJ-1
            KK = KKK
         CASE( 2)
            II = III
            JJ = JJJ+1
            KK = KKK
         CASE(-3)
            II = III
            JJ = JJJ
            KK = KKK-1
         CASE( 3)
            II = III
            JJ = JJJ
            KK = KKK+1
      END SELECT

      ! If the cell is outside the computational domain, check if it is in another mesh

      IF (II<1 .OR. II>M%IBAR .OR. JJ<1 .OR. JJ>M%JBAR .OR. KK<1 .OR. KK>M%KBAR) THEN
         CALL SEARCH_OTHER_MESHES(M%XC(II),M%YC(JJ),M%ZC(KK),NOM,IIO,JJO,KKO)
         IF (NOM>0) M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
         CYCLE SEARCH_LOOP
      ENDIF

      ! Look for thin obstructions bordering the current cell

      DO NN=1,M%N_OBST
         OB=>M%OBSTRUCTION(NN)
         SELECT CASE(IOR)
            CASE(-1)
               IF (II==  OB%I1 .AND. II==  OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE( 1)
               IF (II-1==OB%I1 .AND. II-1==OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE(-2)
               IF (JJ==  OB%J1 .AND. JJ==  OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE( 2)
               IF (JJ-1==OB%J1 .AND. JJ-1==OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE(-3)
               IF (KK==  OB%K1 .AND. KK==  OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
            CASE( 3)
               IF (KK-1==OB%K1 .AND. KK-1==OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
         END SELECT
      ENDDO

      ! If an obstruction is found, assign its cells the current ZONE, just in case the obstruction is removed

      IC = M%CELL_INDEX(II,JJ,KK)
      IF (M%SOLID(IC) .AND. M%OBST_INDEX_C(IC)>0) THEN
         OB => M%OBSTRUCTION(M%OBST_INDEX_C(IC))
         M%PRESSURE_ZONE(OB%I1+1:OB%I2,OB%J1+1:OB%J2,OB%K1+1:OB%K2) = I_ZONE
         CYCLE SEARCH_LOOP
      ENDIF

      ! If the neighboring cell is not solid, assign the pressure zone

      IF (.NOT.M%SOLID(IC) .AND. M%PRESSURE_ZONE(II,JJ,KK)<1) THEN
         Q_N      = Q_N+1
         Q_I(Q_N) = II
         Q_J(Q_N) = JJ
         Q_K(Q_N) = KK
         M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
      ENDIF

   ENDDO SEARCH_LOOP

END DO SORT_QUEUE

DEALLOCATE(Q_I)
DEALLOCATE(Q_J)
DEALLOCATE(Q_K)

END SUBROUTINE ASSIGN_PRESSURE_ZONE



SUBROUTINE GET_N_LAYER_CELLS(DIFFUSIVITY,THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR,N_CELLS,DX_MIN)

! Get number of wall cells in the layer

INTEGER, INTENT(OUT)  :: N_CELLS
REAL(EB), INTENT(OUT) :: DX_MIN
REAL(EB), INTENT(IN)  :: DIFFUSIVITY,THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR
REAL(EB) :: DDSUM
INTEGER  :: N, I

IF (THICKNESS.EQ.0._EB) THEN
   N_CELLS = 0
   DX_MIN = 0._EB
   RETURN
ENDIF
SHRINK_LOOP: DO N=1,999
   DDSUM = 0._EB
   SUM_LOOP: DO I=1,N
      DDSUM = DDSUM + STRETCH_FACTOR**(MIN(I-1,N-I))
   ENDDO SUM_LOOP
   IF ((THICKNESS/DDSUM < CELL_SIZE_FACTOR*SQRT(DIFFUSIVITY)) .OR. (N==999)) THEN
      N_CELLS = N
      DX_MIN = THICKNESS/DDSUM
      EXIT SHRINK_LOOP
   ENDIF
ENDDO SHRINK_LOOP

END SUBROUTINE GET_N_LAYER_CELLS



SUBROUTINE GET_WALL_NODE_COORDINATES(N_CELLS,N_LAYERS,N_LAYER_CELLS,SMALLEST_CELL_SIZE,STRETCH_FACTOR,X_S)

! Get the wall internal coordinates

INTEGER, INTENT(IN)     :: N_CELLS,N_LAYERS, N_LAYER_CELLS(N_LAYERS)
REAL(EB), INTENT(IN)    :: SMALLEST_CELL_SIZE(N_LAYERS),STRETCH_FACTOR
REAL(EB), INTENT(OUT)   :: X_S(0:N_CELLS)

INTEGER I, II, NL
REAL(EB) DX_S

   II = 0
   X_S(0) = 0._EB
   DO NL=1,N_LAYERS
      DO I=1,N_LAYER_CELLS(NL)
         II = II+1
         DX_S = SMALLEST_CELL_SIZE(NL)*STRETCH_FACTOR**(MIN(I-1,N_LAYER_CELLS(NL)-I))
         X_S(II) = X_S(II-1) + DX_S
      ENDDO
   ENDDO

END SUBROUTINE GET_WALL_NODE_COORDINATES



SUBROUTINE GET_WALL_NODE_WEIGHTS(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
         LAYER_THICKNESS,GEOMETRY,X_S,LAYER_DIVIDE,DX,RDX,RDXN,DX_WGT,DXF,DXB,LAYER_INDEX,MF_FRAC)

! Get the wall internal coordinates

INTEGER, INTENT(IN)     :: N_CELLS, N_LAYERS, N_LAYER_CELLS(N_LAYERS), GEOMETRY
REAL(EB), INTENT(IN)    :: X_S(0:N_CELLS),LAYER_THICKNESS(1:N_LAYERS), LAYER_DIVIDE
INTEGER, INTENT(OUT)    :: LAYER_INDEX(0:N_CELLS+1)
REAL(EB), INTENT(OUT)   :: DX(1:N_CELLS),RDX(0:N_CELLS+1),RDXN(0:N_CELLS),DX_WGT(0:N_CELLS),DXF,DXB, &
                           MF_FRAC(1:N_CELLS)

INTEGER I, II, NL, I_GRAD
REAL(EB) R, THICKNESS, X_DIVIDE

   THICKNESS = SUM(LAYER_THICKNESS)

   SELECT CASE(GEOMETRY)
   CASE(SURF_CARTESIAN)
      I_GRAD = 0
   CASE(SURF_CYLINDRICAL)
      I_GRAD = 1
   CASE(SURF_SPHERICAL)
      I_GRAD = 2
   END SELECT

   II = 0
   DO NL=1,N_LAYERS
      DO I=1,N_LAYER_CELLS(NL)
         II = II + 1
         LAYER_INDEX(II) = NL
      ENDDO
   ENDDO
   LAYER_INDEX(0) = 1
   LAYER_INDEX(N_CELLS+1) = N_LAYERS
   DXF = X_S(1)       - X_S(0)
   DXB = X_S(N_CELLS) - X_S(N_CELLS-1)

   ! Compute dx_weight for each node (dx_weight is the ratio of cell size to the combined size of the current and next cell)

   DO I=1,N_CELLS-1
      DX_WGT(I) = (X_S(I)-X_S(I-1))/(X_S(I+1)-X_S(I-1))
   ENDDO
   DX_WGT(0)       = 0.5_EB
   DX_WGT(N_CELLS) = 0.5_EB

   ! Compute dx and 1/dx for each node (dx is the distance from cell boundary to cell boundary)

   DO I=1,N_CELLS
      DX(I)  = X_S(I)-X_S(I-1)
      RDX(I) = 1._EB/DX(I)
   ENDDO

   ! Adjust 1/dx_n to 1/rdr for cylindrical case and 1/r^2dr for spaherical

   IF (GEOMETRY /=SURF_CARTESIAN) THEN
      DO I=1,N_CELLS
         R = THICKNESS-0.5_EB*(X_S(I)+X_S(I-1))
         RDX(I) = RDX(I)/R**I_GRAD
      ENDDO
   ENDIF
   RDX(0)         = RDX(1)
   RDX(N_CELLS+1) = RDX(N_CELLS)

   ! Compute 1/dx_n for each node (dx_n is the distance from cell center to cell center)

   DO I=1,N_CELLS-1
      RDXN(I) = 2._EB/(X_S(I+1)-X_S(I-1))
   ENDDO
   RDXN(0)       = 1._EB/(X_S(1)-X_S(0))
   RDXN(N_CELLS) = 1._EB/(X_S(N_CELLS)-X_S(N_CELLS-1))

   ! Adjust 1/dx_n to r/dr for cylindrical case and r^2/dr for spaherical

   IF (GEOMETRY /= SURF_CARTESIAN) THEN
      DO I=0,N_CELLS
         R = THICKNESS-X_S(I)
         RDXN(I) = RDXN(I)*R**I_GRAD
      ENDDO
   ENDIF

   ! Compute mass flux fraction array (array numbers indicate the fraction of mass flux that is added to the front

   IF (LAYER_DIVIDE >= REAL(N_LAYERS,EB)) THEN
      MF_FRAC = 1.0_EB
   ELSE
      MF_FRAC = 0._EB

      X_DIVIDE = 0._EB
      DO NL=1,N_LAYERS
         IF (LAYER_DIVIDE>=REAL(NL,EB)) THEN
            X_DIVIDE  = X_DIVIDE + LAYER_THICKNESS(NL)
         ELSE
            X_DIVIDE  = X_DIVIDE + MOD(LAYER_DIVIDE,1.0_EB)*LAYER_THICKNESS(NL)
            EXIT
         ENDIF
      ENDDO

      II = 0
      DIVILOOP: DO NL=1,N_LAYERS
         DO I=1,N_LAYER_CELLS(NL)
            II = II + 1
            IF (X_S(II) .LT. X_DIVIDE) THEN
               MF_FRAC(II) = 1._EB
            ELSEIF (X_S(II-1) .LT. X_DIVIDE) THEN
               MF_FRAC(II) = (X_DIVIDE-X_S(II-1))/DX(II)
               EXIT DIVILOOP
            ENDIF
         ENDDO
      ENDDO DIVILOOP
   ENDIF

END SUBROUTINE GET_WALL_NODE_WEIGHTS



SUBROUTINE GET_INTERPOLATION_WEIGHTS(N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS,N_LAYER_CELLS_NEW,X_S,X_S_NEW,INT_WGT)

INTEGER, INTENT(IN)  :: N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS(N_LAYERS),N_LAYER_CELLS_NEW(N_LAYERS)
REAL(EB), INTENT(IN) :: X_S(0:NWP), X_S_NEW(0:NWP_NEW)
REAL(EB), INTENT(OUT) :: INT_WGT(:,:)

REAL(EB) XUP,XLOW,XUP_NEW,XLOW_NEW,DX_NEW
INTEGER I, J, II, JJ, I_BASE, J_BASE, J_OLD,N
II = 0
JJ = 0
I_BASE = 0
J_BASE = 0

INT_WGT = 0._EB
DO N = 1,N_LAYERS
   J_OLD = 1
   DO I = 1,N_LAYER_CELLS_NEW(N)
      II       = I_BASE + I
      XUP_NEW  = X_S_NEW(II)
      XLOW_NEW = X_S_NEW(II-1)
      DX_NEW   = XUP_NEW - XLOW_NEW 
      DO J = J_OLD,N_LAYER_CELLS(N)
         JJ = J_BASE + J
         XUP =  X_S(JJ)
         XLOW = X_S(JJ-1)
         INT_WGT(II,JJ) = (MIN(XUP,XUP_NEW)-MAX(XLOW,XLOW_NEW))/DX_NEW
         IF (XUP .GE. XUP_NEW) EXIT
      ENDDO
      J_OLD = J
   ENDDO
   I_BASE = I_BASE + N_LAYER_CELLS_NEW(N)
   J_BASE = J_BASE + N_LAYER_CELLS(N)
ENDDO

END SUBROUTINE GET_INTERPOLATION_WEIGHTS



SUBROUTINE INTERPOLATE_WALL_ARRAY(NWP,NWP_NEW,INT_WGT,ARR)

INTEGER, INTENT(IN)  :: NWP,NWP_NEW
REAL(EB), INTENT(IN) :: INT_WGT(:,:)
REAL(EB) ARR(NWP),TMP(NWP)

INTEGER I,J

TMP = ARR
ARR = 0._EB
DO I = 1,NWP_NEW
DO J = 1,NWP
   ARR(I) = ARR(I) + INT_WGT(I,J)*TMP(J)
ENDDO
ENDDO

END SUBROUTINE INTERPOLATE_WALL_ARRAY


SUBROUTINE RANDOM_RECTANGLE(XX,YY,ZZ,X1,X2,Y1,Y2,Z1,Z2)

! Choose a random point (XX,YY,ZZ) from within a rectangular volume bounded by X1, X2, ...

REAL(EB), INTENT(IN) :: X1,X2,Y1,Y2,Z1,Z2
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL(EB) :: RN

CALL RANDOM_NUMBER(RN)
XX = X1 + RN*(X2-X1)
CALL RANDOM_NUMBER(RN)
YY = Y1 + RN*(Y2-Y1)
CALL RANDOM_NUMBER(RN)
ZZ = Z1 + RN*(Z2-Z1)

END SUBROUTINE RANDOM_RECTANGLE


SUBROUTINE RANDOM_CONE(XX,YY,ZZ,X0,Y0,Z0,DD,HH)

! Choose a random point (XX,YY,ZZ) from within a vertically oriented cone

REAL(EB), INTENT(IN) :: X0,Y0,Z0,DD,HH
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL(EB) :: RN,THETA,RR

CALL RANDOM_NUMBER(RN)
ZZ = Z0 + HH*RN**2
CALL RANDOM_NUMBER(RN)
THETA = TWOPI*RN
CALL RANDOM_NUMBER(RN)
RR = 0.5_EB*DD*RN*(1._EB-(ZZ-Z0)/HH)
XX = X0 + RR*COS(THETA)
YY = Y0 + RR*SIN(THETA)

END SUBROUTINE RANDOM_CONE


END MODULE GEOMETRY_FUNCTIONS 


 
MODULE PHYSICAL_FUNCTIONS

! Functions for physical quantities

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
IMPLICIT NONE

CONTAINS


SUBROUTINE GET_MASS_FRACTION(Y_IN,INDEX,Y_OUT)

! Y_OUT returns the mass fraction of species INDEX

INTEGER, INTENT(IN) :: INDEX
REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: Y_OUT

Y_OUT = Y2Y_C(INDEX) + DOT_PRODUCT(Y2Y(INDEX,:), Y_IN)
Y_OUT = MIN(1._EB,MAX(0._EB,Y_OUT))

END SUBROUTINE GET_MASS_FRACTION



SUBROUTINE GET_MASS_FRACTION_ALL(Y_IN,Y_OUT)

! Y_OUT returns the mass fraction of all mixture fraction species

REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: Y_OUT(1:N_Y_ARRAY)
INTEGER :: I

DO I=1,N_Y_ARRAY
   Y_OUT(I) = Y2Y_C(I) + DOT_PRODUCT(Y2Y(I,:),Y_IN)
ENDDO

Y_OUT = MIN(1._EB,MAX(0._EB,Y_OUT))

END SUBROUTINE GET_MASS_FRACTION_ALL



SUBROUTINE GET_MOLECULAR_WEIGHT(Y_IN, MW_OUT)

REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: MW_OUT

MW_OUT =  1._EB/ (MW_AVG_Y_C  + DOT_PRODUCT(MW_AVG_Y,Y_IN))

END SUBROUTINE GET_MOLECULAR_WEIGHT



SUBROUTINE GET_SPECIFIC_GAS_CONSTANT(Y_IN, RSUM_OUT)

REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: RSUM_OUT

RSUM_OUT =  R0 * (MW_AVG_Y_C  + DOT_PRODUCT(MW_AVG_Y,Y_IN))

END SUBROUTINE GET_SPECIFIC_GAS_CONSTANT



SUBROUTINE GET_DIFFUSIVITY(Z_IN,YY_SUM,D_OUT,ITMP)

! Compute the diffusivity in the mixture fraction formulation

INTEGER, INTENT(IN) :: ITMP
REAL(EB), INTENT(IN) :: Z_IN(1:I_Z_MAX - I_Z_MIN + 1),YY_SUM
REAL(EB), INTENT(OUT) :: D_OUT
REAL(EB) :: Z(1:I_Z_MAX - I_Z_MIN + 1),OMYYSUM

IF (YY_SUM >=1._EB) THEN
   D_OUT = SPECIES(0)%D(ITMP)
   RETURN
ELSE
   OMYYSUM = 1._EB-MAX(0._EB,YY_SUM)
   Z = MAX(0._EB,MIN(1._EB,Z_IN))/OMYYSUM
ENDIF

D_OUT = (Y2D_C(ITMP) + DOT_PRODUCT(Y2D(ITMP,1:N_MIX_SPECIES),Z))/(Y_MF_SUM_Y_C + DOT_PRODUCT(Y_MF_SUM_Y(I_Z_MIN:I_Z_MAX),Z))

END SUBROUTINE GET_DIFFUSIVITY



SUBROUTINE GET_SPECIFIC_HEAT(Y_IN,CP_OUT,ITMP)

INTEGER, INTENT(IN) :: ITMP
REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: CP_OUT

CP_OUT = Y2CP_C(ITMP) + DOT_PRODUCT(Y2CP(ITMP,:),Y_IN)

END SUBROUTINE GET_SPECIFIC_HEAT



SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT(Y_IN,CPBAR_OUT,ITMP)

INTEGER, INTENT(IN) :: ITMP
REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: CPBAR_OUT

CPBAR_OUT = Y2CPBAR_C(ITMP) + DOT_PRODUCT(Y2CPBAR(ITMP,:),Y_IN)

END SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT



SUBROUTINE GET_SPECIFIC_ENTHALPY(Y_IN,H_G_OUT,ITMP)

INTEGER, INTENT(IN) :: ITMP
REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: H_G_OUT

H_G_OUT = Y2H_G_C(ITMP) + DOT_PRODUCT(Y2H_G(ITMP,:),Y_IN)

END SUBROUTINE GET_SPECIFIC_ENTHALPY



SUBROUTINE GET_CONDUCTIVITY(Y_IN,K_OUT,ITMP)

INTEGER, INTENT(IN) :: ITMP
REAL(EB), INTENT(INOUT) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: K_OUT
REAL(EB) :: MW_MF

CALL GET_MOLECULAR_WEIGHT(Y_IN,MW_MF)
K_OUT = (Y2K_C(ITMP) + DOT_PRODUCT(Y2K(ITMP,:),Y_IN))*MW_MF

END SUBROUTINE GET_CONDUCTIVITY



SUBROUTINE GET_VISCOSITY(Y_IN,MU_OUT,ITMP)

INTEGER, INTENT(IN) :: ITMP
REAL(EB), INTENT(INOUT) :: Y_IN(1:N_SPECIES)
REAL(EB), INTENT(OUT) :: MU_OUT
REAL(EB) :: MW_MF

CALL GET_MOLECULAR_WEIGHT(Y_IN,MW_MF)
MU_OUT = (Y2MU_C(ITMP)+DOT_PRODUCT(Y_IN,Y2MU(ITMP,:)))*MW_MF

END SUBROUTINE GET_VISCOSITY



REAL(EB) FUNCTION DRAG(RE,DRAG_LAW)
 
! drag coefficient

INTEGER, INTENT(IN) :: DRAG_LAW
REAL(EB), INTENT(IN) :: RE

SELECT CASE(DRAG_LAW)

   CASE(SPHERE_DRAG)
      IF (RE<=1._EB) THEN
         DRAG = 24._EB/RE
      ELSEIF (RE>1._EB .AND. RE<1000._EB) THEN
         !!DRAG = 24._EB*(1._EB+0.15_EB*RE**0.687_EB)/RE
         DRAG = 24._EB*(0.85_EB+0.15_EB*RE**0.687_EB)/RE ! smoother curve (RJM)
      ELSEIF (RE>=1000._EB) THEN
         DRAG = 0.44_EB
      ENDIF
      
   CASE(CYLINDER_DRAG)
      IF (RE<=1._EB) THEN
         DRAG = 10._EB/(RE**0.8_EB)
      ELSEIF (RE>1._EB .AND. RE<1000._EB) THEN
         DRAG = 10._EB*(0.6_EB+0.4_EB*RE**0.8_EB)/RE
      ELSEIF (RE>=1000._EB) THEN
         DRAG = 1._EB
      ENDIF
      
END SELECT
 
END FUNCTION DRAG



REAL(EB) FUNCTION SURFACE_DENSITY(NM,IWX,MODE)

! Compute the surface density of a wall cell. 
! If MODE=0, return kg/m2. 
! If MODE=1, return kg/m3. 

INTEGER, INTENT(IN) :: IWX,NM,MODE
INTEGER :: I_GRAD,NWP,II2
REAL(EB) :: THICKNESS,WGT
TYPE (MESH_TYPE), POINTER :: M=>NULL()

M => MESHES(NM)

IF (.NOT.SURFACE(M%IJKW(5,IWX))%THERMALLY_THICK) THEN
   SURFACE_DENSITY = 0._EB
ELSE
   SELECT CASE(SURFACE(M%IJKW(5,IWX))%GEOMETRY)
   CASE(SURF_CARTESIAN)
      I_GRAD = 1
   CASE(SURF_CYLINDRICAL)
      I_GRAD = 2
   CASE(SURF_SPHERICAL)
      I_GRAD = 3
   END SELECT
   IF (SURFACE(M%IJKW(5,IWX))%SHRINK) THEN
      THICKNESS = SUM(M%WALL(IWX)%LAYER_THICKNESS)
      NWP = SUM(M%WALL(IWX)%N_LAYER_CELLS)
      X_S_NEW(0:NWP) = THICKNESS-M%WALL(IWX)%X_S(0:NWP)
   ELSE
      THICKNESS = SUM(SURFACE(M%IJKW(5,IWX))%LAYER_THICKNESS)
      NWP = SURFACE(M%IJKW(5,IWX))%N_CELLS
      X_S_NEW(0:NWP) = THICKNESS-SURFACE(M%IJKW(5,IWX))%X_S(0:NWP)
   ENDIF
   SURFACE_DENSITY = 0._EB
   DO II2=1,NWP
      IF (MODE==0) THEN
         WGT = (X_S_NEW(II2-1)**I_GRAD-X_S_NEW(II2)**I_GRAD)/(REAL(I_GRAD,EB)*SURFACE(M%IJKW(5,IWX))%THICKNESS**(I_GRAD-1))
      ELSE
         WGT = (X_S_NEW(II2-1)**I_GRAD-X_S_NEW(II2)**I_GRAD)/SURFACE(M%IJKW(5,IWX))%THICKNESS**I_GRAD
      ENDIF
      SURFACE_DENSITY = SURFACE_DENSITY + SUM(M%WALL(IWX)%RHO_S(II2,:))*WGT
   ENDDO
ENDIF

END FUNCTION SURFACE_DENSITY



SUBROUTINE DROPLET_SIZE_DISTRIBUTION(DM,RR,CNF,NPT,GAMMA,SIGMA)
 
! Compute droplet Cumulative Number Fraction (CNF)

! CNF(d) = (3/(4*pi))*int_0^d(f(d')*(d'/2)^-3)dd', where f(d') is the PDF of size distribution.
! CNF(di) \approx C * sum_0^i(f(d')*di^-3*ddi)
 
REAL(EB), INTENT(IN) :: DM,GAMMA,SIGMA
INTEGER, INTENT(IN) :: NPT
REAL(EB) :: SUM1,DD1,DI,ETRM,GFAC,SFAC
INTEGER  :: J
REAL(EB), INTENT(OUT) :: RR(0:NPT),CNF(0:NPT)
 
RR(0)  = 0._EB
CNF(0) = 0._EB
SUM1   = 0._EB
! Discretize range of droplet diameters into NPT parts
DD1    = (-LOG(1._EB-0.99_EB)/0.693_EB)**(1._EB/GAMMA)*DM/REAL(NPT,EB)
GFAC   = 0.693_EB*GAMMA*DD1/(DM**GAMMA)
SFAC   = DD1/(SQRT(TWOPI)*SIGMA)
 
INTLOOP: DO J=1,NPT
   DI = (J-.5_EB)*DD1
   RR(J) = .5_EB*DI
   IF (DI<=DM) THEN
      ETRM = EXP(-(LOG(DI/DM))**2/(2._EB*SIGMA**2))
      SUM1 = SUM1 + (SFAC/DI**4)*ETRM
   ELSE
      ETRM = EXP(-0.693_EB*(DI/DM)**GAMMA)
      SUM1 = SUM1 + GFAC*DI**(GAMMA-4._EB)*ETRM
   ENDIF
   CNF(J) = SUM1
ENDDO INTLOOP
 
CNF = CNF/SUM1
 
END SUBROUTINE DROPLET_SIZE_DISTRIBUTION



REAL(EB) FUNCTION FED(Y_IN,RSUM)

! Returns the integrand of FED (Fractional Effective Dose) calculation.

REAL(EB), INTENT(IN) :: Y_IN(1:N_SPECIES),RSUM
REAL(EB) Y_MF_INT, tmp_1

! This is the part depending on gaseous compounds CO2, CO, O2
! Note: Purser uses minutes, here dt is in seconds
!       fed_dose = fed_lco*fed_vco2 + fed_lo

! Next is for CO (ppm)
! CO:  (3.317E-5*RMV*t)/D
!      [RMV]=ltr/min, D=30% COHb concentration at incapacitation
IF (CO_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,CO_INDEX,Y_MF_INT)
   tmp_1 = RCON_MF(CO_INDEX)*Y_MF_INT*1.E6_EB/RSUM
   FED   = 3.317E-5_EB*25.0_EB* tmp_1**(1.036_EB)/(30.0_EB*60.0_EB)
ENDIF
! Next is for CO2
! VCO2: CO2-induced hyperventilation
!      exp(0.1903*c_CO2(%) + 2.0004)
IF (CO2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,CO2_INDEX,Y_MF_INT)
   tmp_1 = RCON_MF(CO2_INDEX)*Y_MF_INT/RSUM
   FED = FED * Exp( 0.1903_EB*tmp_1*100.0_EB + 2.0004_EB )/7.1_EB
ENDIF
! Next is for O2
! LO: low oxygen
! t_Io = exp(8.13-0.54(20.9-%O2)), time in minutes
! F_Io = dt/t_Io
IF (O2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,O2_INDEX,Y_MF_INT)
   tmp_1 = RCON_MF(O2_INDEX)*Y_MF_INT/RSUM
   If ( tmp_1 < 0.20_EB ) FED = FED + 1.0_EB  / (60.0_EB*Exp(8.13_EB-0.54_EB*(20.9_EB-100.0_EB*tmp_1)) )
ENDIF
END FUNCTION FED

SUBROUTINE JANAF_TABLE (I_TMP,CP,H,SPEC_ID,RCON)
!CP data in this subroutine is taken from 4th Edition of the NIST-JANAF Thermochemical Tables
!The tabulated data in the tables have been curve fit
!For gas species that result from liquid evaporation, H is set at the 0 K reference temperature to obtain the correct H_V
!at the H_V reference temperature for that liquid.
!For each species CP is first computed in kJ/mol-K and then converted to kJ/kg-K
CHARACTER(30),INTENT(IN) ::SPEC_ID
INTEGER,INTENT(IN) :: I_TMP
REAL(EB),INTENT(IN) :: RCON
REAL(EB),INTENT(OUT) :: CP,H
REAL(EB) :: CP2,TE

TE = REAL(I_TMP,EB)

SELECT CASE (SPEC_ID)
   CASE DEFAULT
      CP = RCON*GAMMA/(GAMMA-1._EB) !kJ/kg-K
      H = 0._EB !kJ/kg
      
   CASE('AIR')
      IF (TE<100._EB) THEN
         CP = 1032.714_EB + TE * 0.064983_EB !kJ/kg-K
         CP2 = 904.8438_EB + TE * 0.060737_EB  !kJ/kg-K
      ELSE
         CP = 87.05047_EB-18075.59_EB/TE**2 - 669.4452/TE + 0.03857225_EB*TE - 0.00001590889_EB*TE**2 + 0.000000003770398_EB*TE**3 &
              - 0.000000000000459012_EB*TE**4 + 2.232133E-17_EB*TE**5 - 11.54189_EB*LOG(TE) !kJ/mol-K
         CP = CP / MW_N2 * 1000._EB !kJ/kg-K
         CP2 = 72.69494_EB - 566.7158_EB/TE + 0.03684584_EB*TE - 0.0000168704_EB*TE**2 + 0.000000004308058_EB*TE**3 &
              - 0.0000000000005561213_EB*TE**4 + 2.857776E-17_EB*TE**5 - 8.989582_EB*LOG(TE) !kJ/mol-K
         CP2 = CP2 / MW_O2 * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      CP = (1._EB-Y_O2_INFTY)*CP+Y_O2_INFTY*CP2 !kJ/kg-K

   CASE('NITROGEN')
      IF (TE<100._EB) THEN
         CP = 1032.714_EB + TE * 0.064983_EB !kJ/kg-K
      ELSE
         CP = 87.05047_EB-18075.59_EB/TE**2 - 669.4452/TE + 0.03857225_EB*TE - 0.00001590889_EB*TE**2 + 0.000000003770398_EB*TE**3 &
              - 0.000000000000459012_EB*TE**4 + 2.232133E-17_EB*TE**5 - 11.54189_EB*LOG(TE) !kJ/mol-K
         CP = CP / MW_N2 * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      
   CASE('OXYGEN')
      IF (TE<100._EB) THEN
         CP = 904.8438_EB + TE * 0.060737_EB !kJ/kg-K
      ELSE
         CP = 72.69494_EB - 566.7158_EB/TE + 0.03684584_EB*TE - 0.0000168704_EB*TE**2 + 0.000000004308058_EB*TE**3 &
              - 0.0000000000005561213_EB*TE**4 + 2.857776E-17_EB*TE**5 - 8.989582_EB*LOG(TE) !kJ/mol-K
         CP = CP / MW_O2 * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      
   CASE('METHANE')
      IF (TE<100._EB) THEN
         CP = 2078.875_EB - TE * 0.0101122_EB !kJ/kg-K
      ELSE
         CP = -420.3769_EB - 993834.8_EB/TE**2 + 24488.36_EB/TE + 0.03226962_EB*TE - 0.00002813934_EB*TE**2 &
              + 0.000000007651289_EB*TE**3 - 0.0000000000009707785*TE**4+4.772684E-17_EB*TE**5 + 66.26671_EB*LOG(TE) !kJ/mol-K
         CP = CP / 16._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB  !kJ/kg
      
   CASE('ETHYLENE')
      IF (TE<100._EB) THEN
         CP = 1187.5_EB + TE * 0.010002_EB !kJ/kg-K
      ELSE
         CP = -765.804_EB - 964594_EB/TE**2 + 30643.63_EB/TE - 0.07157813_EB*TE + 0.000007792221_EB*TE**2 &
              - 0.0000000003326001_EB*TE**3 - 2.706439E-14_EB*TE**4 + 2.510661E-18_EB*TE**5 + 129.46_EB*LOG(TE) !kJ/mol-K
         CP = CP / 28._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      
   CASE('CARBON DIOXIDE')
      IF (TE<100._EB) THEN
         CP = 658._EB + TE * 0.06981_EB !kJ/kg-K
      ELSE
         CP = -137.8814_EB + 3162.309_EB/TE - 0.01967907_EB*TE + 0.000002677873_EB*TE**2 - 0.0000000002017261_EB*TE**3 &
              + 6.084876E-15_EB*TE**4 + 29.8488_EB*LOG(TE) !kJ/mol-K
         CP = CP / MW_CO2 * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      
   CASE('WATER VAPOR')
      IF (TE<100._EB) THEN
         CP = 1804.5_EB + TE * 0.452714_EB !kJ/kg-K
      ELSE
         CP = 102.2585_EB - 1084.688_EB/TE + 0.04410701_EB*TE - 0.00001181179_EB*TE**2 + 0.000000001601012_EB*TE**3 &
              - 8.408214E-14_EB*TE**4 - 13.55214_EB*LOG(TE) !kJ/mol-K
         CP = CP / MW_H2O * 1000._EB !kJ/kg-K
      ENDIF
      H = 1881295.804_EB !kJ/kg

   CASE('CARBON MONOXIDE')
      IF (TE<100._EB) THEN
         CP = 1033.429_EB + TE * 0.090798_EB !kJ/kg-K
      ELSE
         CP = 43.9033_EB - 69.46386_EB/TE + 0.01806059_EB*TE - 0.000006173999_EB*TE**2 + 0.000000001003309_EB*TE**3 &
              - 6.3073984E-14_EB*TE**4 - 3.423053*LOG(TE) !kJ/mol-K
         CP = CP / MW_CO * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      
   CASE('HYDROGEN')
      IF (TE<100._EB) THEN
         CP = 14685._EB - TE * 9.084204_EB !kJ/kg-K
      ELSE
         CP = 138.9978_EB - 2112.074_EB/TE + 0.06201484_EB*TE - 0.00002746675_EB*TE**2 + 0.00000000767361_EB*TE**3 &
              - 0.000000000001112505_EB*TE**4 + 6.209465E-17_EB*TE**5 -20.90222_EB*LOG(TE) !kJ/mol-K
         CP = CP / MW_H2 * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      
   CASE('ARGON')
         CP = 20.786_EB / 0.040_EB !kJ/kg-K
         H = 0._EB !kJ/kg
               
   CASE('HELIUM')
         CP = 20.786_EB / 0.004_EB !kJ/kg-K
         H = 0._EB !kJ/kg
      
   CASE('SOOT')
!      IF (.NOTE.AL203) TEHEN
         IF (TE<100._EB) THEN
            CP = -39.5_EB + TE * 1.790405_EB !kJ/kg-K
         ELSE
            CP = 138.9978_EB - 2112.074_EB/TE + 0.06201484_EB*TE - 0.00002746675_EB*TE**2 + 0.00000000767361_EB*TE**3 &
               - 0.000000000001112505_EB*TE**4 + 6.209465E-17_EB*TE**5 -20.90222_EB*LOG(TE) !kJ/mol-K
            CP = CP / 12._EB * 1000._EB !kJ/kg-K
         ENDIF
         H = 0._EB !kJ/kg-K
!      ELSE
!         IF (TE<100._EB) THEN
!            CP = -39.5_EB + TE * 0.0214849_EB
!         ELSE
!            CP = 138.9978_EB - 2112.074_EB/TE + 0.06201484_EB*TE - 0.00002746675_EB*TE**2 + 0.00000000767361_EB*TE**3 &
!               - 0.000000000001112505_EB*TE**4 + 6.209465E-17_EB*TE**5 -20.90222_EB*LOG(TE)
!            CP = CP / 102._EB
!         ENDIF
!         H = -89.6177E3_EB
!      ENDIF
   CASE('HYDROGEN BROMIDE')
      IF (TE<100._EB) THEN
         CP = 350.925_EB + TE * 0.084659_EB !kJ/kg-K
      ELSE
         CP = 150.5185+83821.41_EB/TE**2-3380.966_EB/TE+0.05210022_EB*TE-0.00001942684_EB*TE**2+0.000000004370106_EB*TE**3&
              -0.0000000000005138871_EB*TE**4+2.433508E-17_EB*TE**5-21.93194_EB*LOG(TE) !kJ/mol-K
         CP = CP / 81._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
      
   CASE('HYDROGEN CYANIDE')
      IF (TE<100._EB) THEN
         CP = 1070.852_EB + TE * 0.156409_EB !kJ/kg-K
      ELSE
         CP = -57.025353_EB+1556.193_EB/TE+0.001382767_EB*TE-0.000001495245_EB*TE**2+0.0000000001314654_EB*TE**3&
               +1.758821E-14_EB*TE**4-2.419441E-18_EB*TE**5+15.34694_EB*LOG(TE) !kJ/mol-K
         CP = CP / 27._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg

   CASE('HYDROGEN CHLORIDE')
      IF (TE<100._EB) THEN
         CP = 785.111_EB + TE * 0.236287_EB !kJ/kg-K
      ELSE
         CP = 204.4656_EB+183314_EB/TE**2-5890.058_EB/TE+0.06083794_EB*TE-0.00002098845_EB*TE**2+0.000000004486574_EB*TE**3&
              -0.0000000000005086074_EB*TE**4+2.342705E-17_EB*TE**5-30.54396_EB*LOG(TE) !kJ/mol-K
         CP = CP / 36._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg

   CASE('HYDROGEN FLOURIDE')
      IF (TE<100._EB) THEN
         CP = 1370.6_EB + TE * 0.858_EB !kJ/kg-K
      ELSEIF (TE>=100_EB .AND. TE<=200._EB) THEN
         CP = 1456.4_EB !kJ/kg-K
      ELSE
         CP = 204.4656_EB+183314_EB/TE**2-5890.058_EB/TE+0.06083794_EB*TE-0.00002098845_EB*TE**2+0.000000004486574_EB*TE**3&
              -0.0000000000005086074_EB*TE**4+2.342705E-17_EB*TE**5-30.54396_EB*LOG(TE) !kJ/mol-K
         CP = CP / 20._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg

   CASE('NITROGEN OXIDE')
      IF (TE<100._EB) THEN
         CP = 722.2609_EB + TE * 0.010559_EB !kJ/kg-K
      ELSE
         CP = -445.5618_EB-614794.8_EB/TE**2+18374.97_EB/TE-0.08599735_EB*TE+0.00002205842_EB*TE**2-0.0000000039020422_EB*TE**3&
              +0.0000000000003916969_EB*TE**4-1.644367E-17_EB*TE**5+79.24727_EB*LOG(TE) !kJ/mol-K
         CP = CP / 46._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg

   CASE('ACROLEIN')
      IF (TE<50._EB) THEN
         CP = 682.3929_EB + TE * 0.00066_EB !kJ/kg-K
      ELSE
         CP = 332.3828_EB-2602.398_EB/TE+0.6547594_EB*TE-0.0005857_EB*TE**2+0.0000002863821_EB*TE**3 &
              -0.00000000005696995_EB*TE**4-70.66684_EB*LOG(TE) !kJ/mol-K
         CP = CP / 56._EB * 1000._EB !kJ/kg-K
      ENDIF
      H = 0._EB !kJ/kg
   CASE('N-HEPTANE')
      TE = MAX(TE,200._EB)
      CP = -3887.019_EB-2849358._EB/TE**2+130708.2_EB/TE-0.6320126_EB*TE+0.0001695335_EB*TE**2-0.00000003154411_EB*TE**3+&
           0.00000000000323268_EB*TE**4-1.374039E-16_EB*TE**5+670.4689_EB*LOG(TE) !kJ/mol-K
      CP = CP / 100._EB * 1000._EB !kJ/kg-K
      H = 201131._EB  !kJ/kg     
   CASE('METHANOL')
      TE = MAX(TE,200._EB)
      CP = -841.6793_EB-1508086_EB/TE**2+39536.25_EB/TE-0.07191862_EB*TE+0.000006024494_EB*TE**2+0.0000000002922874_EB*TE**3-&
           0.0000000000001171797_EB*TE**4+7.4231E-18_EB*TE**5+138.8215_EB*LOG(TE) !kJ/mol-K
      CP = CP / 32._EB * 1000._EB !kJ/kg-K
      H = 1055464._EB !kJ/kg   
   CASE('ETHANOL')
      TE = MAX(TE,50._EB)
      CP = -445.7989_EB-335159.4_EB/TE**2+15246.42_EB/TE+0.0657487_EB*TE-0.00004985012_EB*TE**2+0.00000001383887_EB*TE**3-&
           0.000000000001815222_EB*TE**4+9.222421E-17_EB*TE**5+78.93689_EB*LOG(TE) !kJ/mol-K
      CP = CP / 46._EB * 1000._EB !kJ/kg-K
      H = 857495.9_EB !kJ/kg       
   CASE('TOLUENE')
      TE = MAX(TE,200._EB)
      CP = -1764.869_EB+730089.6_EB/TE**2+37622.14_EB/TE-0.223475_EB*TE+0.00003599891_EB*TE**2-0.00000000395711_EB*TE**3+&
            0.0000000000002350988_EB*TE**4-5.2885E-18_EB*TE**5+315.3972_EB*LOG(TE) !kJ/mol-K
      CP = CP / 92._EB * 1000._EB !kJ/kg-K
      H = 356078.2_EB !kJ/kg      
END SELECT    

END SUBROUTINE JANAF_TABLE

SUBROUTINE JANAF_TABLE_LIQUID (I_TMP,CP,H_V,H_L,T_REF,T_MELT,T_BOIL,SPEC_ID)
!CP data in this subroutine is taken from 4th Edition of the NIST-JANAF Thermochemical Tables
!The tabulated data in the tables have been curve fit
!For each species CP is first computed in kJ/mol-K and then converted to kJ/kg-K
CHARACTER(30), INTENT(IN) ::SPEC_ID
INTEGER, INTENT(IN) :: I_TMP
REAL(EB), INTENT(INOUT) :: CP,H_V,H_L,T_REF,T_MELT,T_BOIL
REAL(EB) :: T

SELECT CASE (SPEC_ID)
   CASE('WATER VAPOR')
      T = MAX(273._EB,REAL(I_TMP,EB))
      CP = 194.74308_EB - 1.2338708_EB*T + 0.0047935684_EB*T**2 - 0.0000084183512_EB*T**3 + 0.0000000057563403_EB*T**4 !kJ/mol-K
      CP = CP / MW_H2O * 1000._EB !kJ/kg-K
      H_L = -1150345._EB!kJ/kg
      H_V = 2385253.83_EB !kJ/kg
      T_REF = 273.15_EB !K
      T_MELT = 273.15_EB !K
      T_BOIL = 373.15_EB !K
   CASE('N-HEPTANE')
      T = MAX(183._EB,REAL(I_TMP,EB))
      CP =  0.001061074_EB*T**2 - 0.267961_EB*T + 209.8049_EB !kJ/mol-K
      CP = CP / 100._EB * 1000._EB !kJ/kg-K
      H_L = -358349.9_EB !kJ/kg
      H_V = 434962.3_EB !kJ/kg
      T_REF = 182.55_EB !K
      T_MELT = 182.55_EB !K
      T_BOIL = 371.53_EB !K
   CASE('METHANOL')
      T = MAX(177._EB,REAL(I_TMP,EB))
      CP =   0.0008750616_EB*T**2 - 0.3274918_EB*T + 101.2359_EB !kJ/mol-K
      CP = CP / 32._EB * 1000._EB !kJ/kg-K
      H_L = -390974.5_EB !kJ/kg-K
      H_V = 1274941.64_EB !kJ/kg
      T_REF = 176.6_EB !K
      T_MELT = 176.6_EB !K
      T_BOIL = 337.7_EB !K
   CASE('ETHANOL')
      T = MAX(159._EB,REAL(I_TMP,EB))
      CP =   0.3229_EB*T+16.272_EB
      CP = CP / 46._EB * 1000._EB !kJ/kg-K
      H_L = -233706.2_EB !kJ/kg-K
      H_V = 1000197._EB !kJ/kg
      T_REF = 158.8_EB !K
      T_MELT = 158.8_EB !K
      T_BOIL = 358.8_EB !K
   CASE('TOLUENE')
      T = MAX(178._EB,REAL(I_TMP,EB))
      CP = 0.0007691129_EB*T**2- 0.1856684_EB*T + 143.1817_EB !kJ/mol-K
      CP = CP / 92._EB * 1000._EB !kJ/kg-K
      H_L = -260230.8_EB !kJ/kg-K
      H_V = 490124._EB !kJ/kg
      T_REF = 178.1_EB !K
      T_MELT = 178.1_EB !K
      T_BOIL = 383.8_EB !K      
   CASE DEFAULT
      CP = -1._EB
      H_L = -1._EB
      H_V = -1._EB
      T_REF = -1._EB
END SELECT    

END SUBROUTINE JANAF_TABLE_LIQUID


SUBROUTINE GAS_PROPS(GAS_NAME,SIGMA,EPSOK,MW,ABSORBING,FORMULA)
 
! Molecular weight (g/mol) and Lennard-Jones properties
 
REAL(EB) SIGMA,EPSOK,MW,SIGMAIN,EPSOKIN,MWIN
CHARACTER(30) GAS_NAME,FORMULA
LOGICAL ABSORBING
 
SIGMAIN = SIGMA
EPSOKIN = EPSOK
MWIN    = MW
 
SELECT CASE(GAS_NAME)
   CASE('AIR')             
      SIGMA=3.711_EB 
      EPSOK= 78.6_EB  
      MW=MW_AIR
      FORMULA='Air'
   CASE('CARBON MONOXIDE') 
      SIGMA=3.690_EB 
      EPSOK= 91.7_EB  
      MW=28._EB
      ABSORBING = .TRUE.
      FORMULA='CO'
   CASE('CARBON DIOXIDE')  
      SIGMA=3.941_EB 
      EPSOK=195.2_EB  
      MW=44._EB
      ABSORBING = .TRUE.      
      FORMULA='CO2'
   CASE('ETHYLENE')        
      SIGMA=4.163_EB 
      EPSOK=224.7_EB 
      MW=28._EB
      ABSORBING = .TRUE.
      FORMULA='C2H2'
   CASE('ARGON')
      SIGMA=3.42_EB
      EPSOK= 124.0_EB
      MW= 40._EB
      FORMULA='Ar'
   CASE('HELIUM')          
      SIGMA=2.551_EB 
      EPSOK= 10.22_EB 
      MW= 4._EB
      FORMULA='He2'
   CASE('HYDROGEN')        
      SIGMA=2.827_EB 
      EPSOK= 59.7_EB
      MW= 2._EB
      FORMULA='H2'
   CASE('METHANE')         
      SIGMA=3.758_EB 
      EPSOK=148.6_EB  
      MW=16._EB
      ABSORBING = .TRUE.
      FORMULA='CH4'
   CASE('NITROGEN')        
      SIGMA=3.798_EB 
      EPSOK= 71.4_EB  
      MW=28._EB
      FORMULA='N2'
   CASE('OXYGEN')          
      SIGMA=3.467_EB 
      EPSOK=106.7_EB  
      MW=32._EB
      FORMULA='O2'
   CASE('PROPANE')         
      SIGMA=5.118_EB 
      EPSOK=237.1_EB  
      MW=44._EB
      ABSORBING = .TRUE.
      FORMULA='C3H8'
   CASE('WATER VAPOR')     
      SIGMA=2.641_EB 
      EPSOK=809.1_EB  
      MW=18._EB
      ABSORBING = .TRUE.
      FORMULA='H2O'
   CASE DEFAULT            
      SIGMA=3.711_EB 
      EPSOK= 78.6_EB
      MW=MW_AIR
   CASE('SOOT')        
      SIGMA=3.798_EB 
      EPSOK= 71.4_EB  
      MW=12._EB
      FORMULA='C'
   CASE('N-HEXANE')
      SIGMA=4.524_EB
      EPSOK=199.41_EB
      MW=86._EB
      FORMULA='C6H14'
   CASE('N-HEPTANE')
      SIGMA=4.701_EB
      EPSOK=205.78_EB
      MW=100._EB
      FORMULA='C7H16'
   CASE('N-OCTANE')
      SIGMA=4.892_EB
      EPSOK=231.16_EB
      MW=114._EB
      FORMULA='C8H18'
   CASE('N-DECANE')
      SIGMA=5.233_EB
      EPSOK=226.46_EB
      MW=142._EB
      FORMULA='C10H22'  
   CASE('METHANOL')
      SIGMA=3.626_EB
      EPSOK=481.8_EB
      MW=32._EB
      FORMULA='CH3OH'  
   CASE('ETHANOL')
      SIGMA=4.530_EB
      EPSOK=362.6_EB
      MW=46._EB
      FORMULA='C2H5OH'  
    CASE('HYDROGEN BROMIDE')   
      SIGMA=3.353_EB
      EPSOK=449._EB
      MW=81._EB
      FORMULA='HBR'  
    CASE('HYDROGEN CHLORIDE')   
      SIGMA=3.339_EB
      EPSOK=344.7_EB
      MW=36._EB
      FORMULA='HCL'  
    CASE('HYDROGEN FLOURIDE')   
      SIGMA=3.148_EB
      EPSOK=330._EB
      MW=20._EB
      FORMULA='HF'  
    CASE('HYDROGEN CYANIDE')   
      SIGMA=3.63_EB
      EPSOK=569.1_EB
      MW=27._EB
      FORMULA='HCN'  
    CASE('TOLUENE')   
      SIGMA=5.698_EB
      EPSOK=480._EB !Estimated from BENZENE
      MW=92._EB
      FORMULA='HCN'  
END SELECT  
 
IF (SIGMAIN>0._EB) SIGMA = SIGMAIN
IF (EPSOKIN>0._EB) EPSOK = EPSOKIN
IF (MWIN   >0._EB) MW    = MWIN 
 
END SUBROUTINE GAS_PROPS


END MODULE PHYSICAL_FUNCTIONS

 
MODULE MATH_FUNCTIONS

USE PRECISION_PARAMETERS
IMPLICIT NONE 
 
CONTAINS
 
REAL(EB) FUNCTION AFILL(A111,A211,A121,A221,A112,A212,A122,A222,P,R,S)
! Linear interpolation function
REAL(EB) :: A111,A211,A121,A221,A112,A212,A122,A222,P,R,S,PP,RR,SS
PP = 1._EB-P
RR = 1._EB-R
SS = 1._EB-S
AFILL = ((PP*A111+P*A211)*RR+(PP*A121+P*A221)*R)*SS+((PP*A112+P*A212)*RR+(PP*A122+P*A222)*R)*S
END FUNCTION AFILL
 

REAL(EB) FUNCTION AFILL2(A,I,J,K,P,R,S)
! Linear interpolation function. Same as AFILL, only it reads in entire array.
REAL(EB), INTENT(IN), DIMENSION(0:,0:,0:) :: A
INTEGER, INTENT(IN) :: I,J,K
REAL(EB) A111,A211,A121,A221,A112,A212,A122,A222,P,R,S,PP,RR,SS
A111 = A(I,J,K)
A211 = A(I+1,J,K)
A121 = A(I,J+1,K)
A221 = A(I+1,J+1,K)
A112 = A(I,J,K+1)
A212 = A(I+1,J,K+1)
A122 = A(I,J+1,K+1)
A222 = A(I+1,J+1,K+1)
PP = 1._EB-P
RR = 1._EB-R
SS = 1._EB-S
AFILL2 = ((PP*A111+P*A211)*RR+(PP*A121+P*A221)*R)*SS+ ((PP*A112+P*A212)*RR+(PP*A122+P*A222)*R)*S
END FUNCTION AFILL2


REAL(EB) FUNCTION POLYVAL(N,TEMP,COEF)
! Calculate the value of polynomial function.
INTEGER N,I
REAL(EB) TEMP, COEF(N), VAL
VAL = 0._EB
DO I=1,N
   VAL  = VAL  + COEF(I)*TEMP**(I-1)
ENDDO
POLYVAL = VAL
END FUNCTION POLYVAL
 
 
SUBROUTINE GET_RAMP_INDEX(ID,TYPE,RAMP_INDEX)
USE GLOBAL_CONSTANTS, ONLY: N_RAMP,RAMP_ID,RAMP_TYPE
CHARACTER(*), INTENT(IN) :: ID,TYPE
INTEGER, INTENT(OUT) :: RAMP_INDEX
INTEGER :: NR
 
IF (ID=='null') THEN
   RAMP_INDEX = 0
   RETURN
ENDIF
 
SEARCH: DO NR=1,N_RAMP
   IF (ID==RAMP_ID(NR)) THEN
      RAMP_INDEX = NR
      RETURN
   ENDIF
ENDDO SEARCH
 
N_RAMP                = N_RAMP + 1
RAMP_INDEX            = N_RAMP
RAMP_ID(RAMP_INDEX)   = ID
RAMP_TYPE(RAMP_INDEX) = TYPE
END SUBROUTINE GET_RAMP_INDEX

SUBROUTINE GET_TABLE_INDEX(ID,TYPE,TABLE_INDEX)
USE GLOBAL_CONSTANTS, ONLY: N_TABLE,TABLE_ID,TABLE_TYPE
CHARACTER(*), INTENT(IN) :: ID
INTEGER, INTENT(IN) :: TYPE
INTEGER, INTENT(OUT) :: TABLE_INDEX
INTEGER :: NT
 
IF (ID=='null') THEN
   TABLE_INDEX = 0
   RETURN
ENDIF
 
SEARCH: DO NT=1,N_TABLE
   IF (ID==TABLE_ID(NT)) THEN
      TABLE_INDEX = NT
      RETURN
   ENDIF
ENDDO SEARCH
 
N_TABLE                = N_TABLE + 1
TABLE_INDEX            = N_TABLE
TABLE_ID(TABLE_INDEX)   = ID
TABLE_TYPE(TABLE_INDEX) = TYPE

END SUBROUTINE GET_TABLE_INDEX



REAL(EB) FUNCTION EVALUATE_RAMP(RAMP_INPUT,TAU,RAMP_INDEX)

! General time ramp up
 
USE TYPES, ONLY: RAMPS
USE DEVICE_VARIABLES, ONLY: DEVICE
USE CONTROL_VARIABLES, ONLY: CONTROL
REAL(EB), INTENT(IN) :: RAMP_INPUT,TAU
REAL(EB):: RAMP_POSITION,VALUE_1,VALUE_2,POSITION_VALUE
INTEGER :: POSITION_INDEX
INTEGER,INTENT(IN)   :: RAMP_INDEX

SELECT CASE(RAMP_INDEX)
   CASE(-2)
      EVALUATE_RAMP = MAX(TANH(RAMP_INPUT/TAU),0._EB)
   CASE(-1)
      EVALUATE_RAMP = MIN( (RAMP_INPUT/TAU)**2 , 1.0_EB )
   CASE( 0)
      EVALUATE_RAMP = 1._EB
   CASE(1:)
      IF (RAMPS(RAMP_INDEX)%DEVC_INDEX > 0) THEN
         RAMP_POSITION = &
            MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,DEVICE(RAMPS(RAMP_INDEX)%DEVC_INDEX)%INSTANT_VALUE - RAMPS(RAMP_INDEX)%T_MIN))
         EVALUATE_RAMP = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(NINT(RAMP_POSITION/RAMPS(RAMP_INDEX)%DT))
      ELSE
         RAMP_POSITION = MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,RAMP_INPUT - RAMPS(RAMP_INDEX)%T_MIN))
         POSITION_VALUE = RAMP_POSITION/RAMPS(RAMP_INDEX)%DT
         POSITION_INDEX = FLOOR(POSITION_VALUE)
         VALUE_1 = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(POSITION_INDEX)
         VALUE_2 = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(POSITION_INDEX+1)
         EVALUATE_RAMP = VALUE_1 + (POSITION_VALUE-POSITION_INDEX)*(VALUE_2-VALUE_1)
      ENDIF
END SELECT
 
END FUNCTION EVALUATE_RAMP

 
REAL(EB) FUNCTION ERFC(X)
 
! Complimentary ERF function
 
REAL(EB), INTENT(IN) :: X
REAL(EB) ERFCS(13), ERFCCS(24), ERC2CS(23),XSML,XMAX,SQEPS,Y
DATA ERFCS( 1) /   -.049046121234691808_EB /
DATA ERFCS( 2) /   -.14226120510371364_EB /
DATA ERFCS( 3) /    .010035582187599796_EB /
DATA ERFCS( 4) /   -.000576876469976748_EB /
DATA ERFCS( 5) /    .000027419931252196_EB /
DATA ERFCS( 6) /   -.000001104317550734_EB /
DATA ERFCS( 7) /    .000000038488755420_EB /
DATA ERFCS( 8) /   -.000000001180858253_EB /
DATA ERFCS( 9) /    .000000000032334215_EB /
DATA ERFCS(10) /   -.000000000000799101_EB /
DATA ERFCS(11) /    .000000000000017990_EB /
DATA ERFCS(12) /   -.000000000000000371_EB /
DATA ERFCS(13) /    .000000000000000007_EB /
DATA ERC2CS( 1) /   -.069601346602309501_EB /
DATA ERC2CS( 2) /   -.041101339362620893_EB /
DATA ERC2CS( 3) /    .003914495866689626_EB /
DATA ERC2CS( 4) /   -.000490639565054897_EB /
DATA ERC2CS( 5) /    .000071574790013770_EB /
DATA ERC2CS( 6) /   -.000011530716341312_EB /
DATA ERC2CS( 7) /    .000001994670590201_EB /
DATA ERC2CS( 8) /   -.000000364266647159_EB /
DATA ERC2CS( 9) /    .000000069443726100_EB /
DATA ERC2CS(10) /   -.000000013712209021_EB /
DATA ERC2CS(11) /    .000000002788389661_EB /
DATA ERC2CS(12) /   -.000000000581416472_EB /
DATA ERC2CS(13) /    .000000000123892049_EB /
DATA ERC2CS(14) /   -.000000000026906391_EB /
DATA ERC2CS(15) /    .000000000005942614_EB /
DATA ERC2CS(16) /   -.000000000001332386_EB /
DATA ERC2CS(17) /    .000000000000302804_EB /
DATA ERC2CS(18) /   -.000000000000069666_EB /
DATA ERC2CS(19) /    .000000000000016208_EB /
DATA ERC2CS(20) /   -.000000000000003809_EB /
DATA ERC2CS(21) /    .000000000000000904_EB /
DATA ERC2CS(22) /   -.000000000000000216_EB /
DATA ERC2CS(23) /    .000000000000000052_EB /
DATA ERFCCS( 1) /     0.0715179310202925_EB /
DATA ERFCCS( 2) /   -.026532434337606719_EB /
DATA ERFCCS( 3) /    .001711153977920853_EB /
DATA ERFCCS( 4) /   -.000163751663458512_EB /
DATA ERFCCS( 5) /    .000019871293500549_EB /
DATA ERFCCS( 6) /   -.000002843712412769_EB /
DATA ERFCCS( 7) /    .000000460616130901_EB /
DATA ERFCCS( 8) /   -.000000082277530261_EB /
DATA ERFCCS( 9) /    .000000015921418724_EB /
DATA ERFCCS(10) /   -.000000003295071356_EB /
DATA ERFCCS(11) /    .000000000722343973_EB /
DATA ERFCCS(12) /   -.000000000166485584_EB /
DATA ERFCCS(13) /    .000000000040103931_EB /
DATA ERFCCS(14) /   -.000000000010048164_EB /
DATA ERFCCS(15) /    .000000000002608272_EB /
DATA ERFCCS(16) /   -.000000000000699105_EB /
DATA ERFCCS(17) /    .000000000000192946_EB /
DATA ERFCCS(18) /   -.000000000000054704_EB /
DATA ERFCCS(19) /    .000000000000015901_EB /
DATA ERFCCS(20) /   -.000000000000004729_EB /
DATA ERFCCS(21) /    .000000000000001432_EB /
DATA ERFCCS(22) /   -.000000000000000439_EB /
DATA ERFCCS(23) /    .000000000000000138_EB /
DATA ERFCCS(24) /   -.000000000000000048_EB /
 
XSML = -200._EB
XMAX = 200._EB
SQEPS = 0.001_EB
 
IF (X<=XSML) THEN
   ERFC = 2._EB
ELSE
 
   IF (X<=XMAX) THEN
      Y = ABS(X)
      IF (Y<=1.0_EB) THEN  ! ERFC(X) = 1.0 - ERF(X) FOR -1._EB <= X <= 1.
         IF (Y<SQEPS)  ERFC = 1.0_EB - 2.0_EB*X/SQRTPI
         IF (Y>=SQEPS) ERFC = 1.0_EB - X*(1.0_EB + CSEVL (2._EB*X*X-1._EB, ERFCS, 10) )
      ELSE  ! ERFC(X) = 1.0 - ERF(X) FOR 1._EB < ABS(X) <= XMAX
         Y = Y*Y
         IF (Y<=4._EB) ERFC = EXP(-Y)/ABS(X) * (0.5_EB + CSEVL ((8._EB/Y-5._EB)/3._EB,ERC2CS, 10) )
         IF (Y>4._EB) ERFC = EXP(-Y)/ABS(X) * (0.5_EB + CSEVL (8._EB/Y-1._EB,ERFCCS, 10) )
         IF (X<0._EB) ERFC = 2.0_EB - ERFC
      ENDIF
   ELSE
      ERFC = 0._EB
   ENDIF
ENDIF
RETURN
 
END FUNCTION ERFC
 
 
SUBROUTINE GAUSSJ(A,N,NP,B,M,MP,IERROR)
 
! Solve a linear system of equations with Gauss-Jordon elimination
! Source: Press et al. "Numerical Recipes"
 
INTEGER :: M,MP,N,NP,I,ICOL,IROW,J,K,L,LL,INDXC(NP),INDXR(NP),IPIV(NP)
REAL(EB) :: A(NP,NP),B(NP,MP),BIG,DUM,PIVINV
INTEGER, INTENT(OUT) :: IERROR
 
IERROR = 0
IPIV(1:N) = 0
 
DO I=1,N
   BIG = 0._EB
   DO J=1,N
      IF (IPIV(J)/=1) THEN
         DO K=1,N
            IF (IPIV(K)==0) THEN
               IF (ABS(A(J,K))>=BIG) THEN
                  BIG = ABS(A(J,K))
                  IROW = J
                  ICOL = K
               ENDIF
            ELSE IF (IPIV(K)>1) THEN
               IERROR = 103   ! Singular matrix in gaussj
               RETURN
            ENDIF
         ENDDO
      ENDIF
   ENDDO
   IPIV(ICOL) = IPIV(ICOL) + 1
   IF (IROW/=ICOL) THEN
      DO L=1,N
         DUM = A(IROW,L)
         A(IROW,L) = A(ICOL,L)
         A(ICOL,L) = DUM
      ENDDO
      DO L=1,M
         DUM = B(IROW,L)
         B(IROW,L) = B(ICOL,L)
         B(ICOL,L) = DUM
      ENDDO
   ENDIF
   INDXR(I) = IROW
   INDXC(I) = ICOL
   IF (A(ICOL,ICOL)==0._EB) THEN
      IERROR = 103  ! Singular matrix in gaussj
      RETURN
      ENDIF
   PIVINV = 1._EB/A(ICOL,ICOL)
   A(ICOL,ICOL) = 1._EB
   A(ICOL,1:N) = A(ICOL,1:N) * PIVINV
   B(ICOL,1:M) = B(ICOL,1:M) * PIVINV
   DO LL=1,N
      IF (LL/=ICOL) THEN
         DUM = A(LL,ICOL)
         A(LL,ICOL) = 0._EB
         A(LL,1:N) = A(LL,1:N) - A(ICOL,1:N)*DUM
         B(LL,1:M) = B(LL,1:M) - B(ICOL,1:M)*DUM
      ENDIF
   ENDDO
ENDDO
DO L=N,1,-1
   IF (INDXR(L)/=INDXC(L)) THEN
      DO K=1,N
         DUM = A(K,INDXR(L))
         A(K,INDXR(L)) = A(K,INDXC(L))
         A(K,INDXC(L)) = DUM
      ENDDO
   ENDIF
ENDDO
 
END SUBROUTINE GAUSSJ
 
 
REAL(EB) FUNCTION CSEVL(X,CS,N)
 
REAL(EB), INTENT(IN) :: X
REAL(EB) CS(:),B1,B0,TWOX,B2
INTEGER NI,N,I
 
B1=0._EB
B0=0._EB
TWOX=2._EB*X
DO I=1,N
B2=B1
B1=B0
NI=N+1-I
B0=TWOX*B1-B2+CS(NI)
ENDDO
 
CSEVL = 0.5_EB*(B0-B2)
 
END FUNCTION CSEVL

INTEGER(2) FUNCTION TWO_BYTE_REAL(REAL_IN)
REAL(FB),INTENT(IN) :: REAL_IN
INTEGER(2) EXP,TEMP

IF (ABS(REAL_IN) <= 1.E-16_8) THEN
   TWO_BYTE_REAL = 0
ELSEIF (ABS(REAL_IN) >= 1.E+16_8) THEN
   TEMP = 999
   EXP=30
   TWO_BYTE_REAL = TEMP * 2**5
   TWO_BYTE_REAL = TWO_BYTE_REAL + EXP
ELSE
   EXP = FLOOR(LOG10(ABS(REAL_IN)))
   TEMP = ABS(REAL_IN * 10**(-REAL(EXP+1,8)+3))
   EXP = EXP + 15
   TWO_BYTE_REAL = TEMP * 2**5
   TWO_BYTE_REAL = TWO_BYTE_REAL + EXP
ENDIF
IF (REAL_IN < 0._8) TWO_BYTE_REAL = IBSET(TWO_BYTE_REAL,15)

END FUNCTION TWO_BYTE_REAL

INTEGER FUNCTION RLE_COMPRESSION(QIN,NQIN,QMIN,QMAX,COUT)

! Compress the array QIN(1:NQIN) by first mapping to one byte integers and then
! using Run-Length Encoding to convert runs of common integers IIIIIIII TO  #In
! where #=255 is a marker character to distinguish between literal characters 
! and runs of n repeats

INTEGER, INTENT(IN) :: NQIN
REAL, INTENT(IN), DIMENSION(NQIN) :: QIN
REAL, INTENT(IN) :: QMIN, QMAX
CHARACTER(1), DIMENSION(NQIN) :: COUT
CHARACTER(1), PARAMETER :: MARK=CHAR(255)
CHARACTER(1) :: THISCHAR, LASTCHAR
INTEGER :: IIN,IOUT,NREPEATS,IQVAL

IF (QMAX<=QMIN .OR. NQIN<=0) THEN
   RLE_COMPRESSION = 0
   RETURN
ENDIF

IOUT=1
LASTCHAR=MARK

DO IIN = 1, NQIN
   IQVAL = 254*(QIN(IIN) - QMIN)/(QMAX-QMIN)
   IF (IQVAL<0)   IQVAL=0
   IF (IQVAL>254) IQVAL=254
   THISCHAR =CHAR(IQVAL)

   IF (THISCHAR == LASTCHAR) THEN
      NREPEATS = NREPEATS + 1
   ELSE
      NREPEATS = 1
   ENDIF

   IF (NREPEATS>=1 .AND. NREPEATS<=3) THEN
      COUT(IOUT) = THISCHAR
      LASTCHAR=THISCHAR
   ELSEIF (NREPEATS>=4) THEN
      IF (NREPEATS==4) THEN
         IOUT = IOUT - 3
         COUT(IOUT) = MARK
         IOUT = IOUT + 1
         COUT(IOUT) = THISCHAR
         IOUT = IOUT + 1
      ENDIF
      IF (NREPEATS/=4) IOUT = IOUT - 1
      COUT(IOUT) = CHAR(NREPEATS)
      IF (NREPEATS==254) THEN
          NREPEATS=1
          LASTCHAR=MARK
      ENDIF
   ENDIF
   IOUT = IOUT + 1
ENDDO

RLE_COMPRESSION = IOUT - 1

END FUNCTION RLE_COMPRESSION

SUBROUTINE INTERPOLATE1D(X,Y,XI,ANS)
REAL(EB), INTENT(IN), DIMENSION(:) :: X, Y
REAL(EB), INTENT(IN) :: XI
REAL(EB), INTENT(OUT) :: ANS
INTEGER I, UX,LX

UX = UBOUND(X,1)
LX = LBOUND(X,1)

IF (XI <= X(LX)) THEN
  ANS = Y(LX)
ELSEIF (XI >= X(UX)) THEN
  ANS = Y(UX)
ELSE
  L1: DO I=LX,UX-1
    IF (XI -X(I) == 0._EB) THEN
      ANS = Y(I)
      EXIT L1
    ELSEIF (X(I+1)>XI) THEN
      ANS = Y(I)+(XI-X(I))/(X(I+1)-X(I)) * (Y(I+1)-Y(I))
      EXIT L1
    ENDIF
  ENDDO L1
ENDIF

END SUBROUTINE INTERPOLATE1D

REAL(EB) FUNCTION NORMAL(MEAN,SIGMA) !returns a normal distribution
IMPLICIT NONE
REAL(EB) :: MEAN,SIGMA,TMP,FAC,GSAVE,RSQ,R1,R2,RN
INTEGER :: FLAG
SAVE FLAG,GSAVE
DATA FLAG /0/
IF (FLAG.EQ.0) THEN
   RSQ=2.0_EB
   DO WHILE(RSQ.GE.1.0_EB.OR.RSQ.EQ.0.0_EB)
      CALL RANDOM_NUMBER(RN)
      R1=2.0_EB*RN-1.0_EB
      CALL RANDOM_NUMBER(RN)
      R2=2.0_EB*RN-1.0_EB
      RSQ=R1*R1+R2*R2
   ENDDO
   FAC=SQRT(-2.0_EB*LOG(RSQ)/RSQ)
   GSAVE=R1*FAC
   TMP=R2*FAC
   FLAG=1
ELSE
   TMP=GSAVE
   FLAG=0
ENDIF
NORMAL=TMP*SIGMA+MEAN
RETURN
END FUNCTION NORMAL

SUBROUTINE CROSS_PRODUCT(C,A,B)
IMPLICIT NONE
REAL(EB), INTENT(IN) :: A(3),B(3)
REAL(EB), INTENT(OUT) :: C(3)
! C = A x B
C(1) = A(2)*B(3)-A(3)*B(2)
C(2) = A(3)*B(1)-A(1)*B(3)
C(3) = A(1)*B(2)-A(2)*B(1)
END SUBROUTINE CROSS_PRODUCT

REAL(EB) FUNCTION NORM2(V)
IMPLICIT NONE
REAL(EB), INTENT(IN) :: V(3)
! Returns the L2 norm of the 3-vector V
NORM2 = SQRT(DOT_PRODUCT(V,V))
RETURN
END FUNCTION


END MODULE MATH_FUNCTIONS


MODULE TRAN 
 
! Coordinate transformation functions
 
USE PRECISION_PARAMETERS
IMPLICIT NONE
TYPE TRAN_TYPE
   REAL(EB), POINTER, DIMENSION(:,:) :: C1=>NULL(),C2=>NULL(),C3=>NULL(),CCSTORE=>NULL(),PCSTORE=>NULL()
   INTEGER, POINTER, DIMENSION(:,:) :: IDERIVSTORE=>NULL()
   INTEGER NOC(3),ITRAN(3),NOCMAX
END TYPE TRAN_TYPE
TYPE (TRAN_TYPE), ALLOCATABLE, TARGET, DIMENSION(:) :: TRANS
 
 
CONTAINS
 
 
REAL(EB) FUNCTION G(X,IC,NM)
 
! Coordinate transformation function
 
REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T=>NULL()
 
T => TRANS(NM)
 
N = T%NOC(IC)
IF (N==0) THEN
   G = X
   RETURN
ENDIF
 
SELECT CASE(T%ITRAN(IC))
   CASE(1)
      G = 0._EB
      DO I=1,N+1
         G = G + T%C1(I,IC)*X**I
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      G = T%C2(II-1,IC) + T%C3(II,IC)*(X-T%C1(II-1,IC))
END SELECT
 
END FUNCTION G
 
 
REAL(EB) FUNCTION GP(X,IC,NM)
 
! Derivative of the coordinate transformation function
 
REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T=>NULL()
 
T => TRANS(NM)
N =  T%NOC(IC)
IF (N==0) THEN
   GP = 1._EB
   RETURN
ENDIF
 
SELECT CASE(T%ITRAN(IC)) 
   CASE(1)
      GP = 0._EB
      DO I=1,N+1
         GP = GP + I*T%C1(I,IC)*X**(I-1)
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GP = T%C3(II,IC)
END SELECT
 
END FUNCTION GP
 
 
REAL(EB) FUNCTION GINV(Z,IC,NM)
 
! Inverse of the coordinate transformation function
 
REAL(EB) :: GF
INTEGER :: N,IT,II,I
REAL(EB), INTENT(IN) :: Z
INTEGER, INTENT(IN)  :: IC,NM
TYPE (TRAN_TYPE), POINTER :: T=>NULL()
 
T => TRANS(NM)
GINV = Z
N = T%NOC(IC)
IF (N==0) RETURN
 
SELECT CASE(T%ITRAN(IC))
   CASE(1)
      LOOP1: DO IT=1,10
         GF = G(GINV,IC,NM)-Z
         IF (ABS(GF)<0.002_EB) EXIT LOOP1
         GINV = GINV - GF/GP(GINV,IC,NM)
      ENDDO LOOP1
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (Z<=T%C2(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GINV = T%C1(II-1,IC) + (Z-T%C2(II-1,IC))/T%C3(II,IC)
END SELECT
 
END FUNCTION GINV
 
 
END MODULE TRAN

MODULE OPENMP

! Module for OpenMP check

USE GLOBAL_CONSTANTS

IMPLICIT NONE

PUBLIC OPENMP_CHECK
 
CONTAINS

SUBROUTINE OPENMP_CHECK

!$ INTEGER :: OMP_GET_NUM_THREADS

!$OMP PARALLEL
!$OMP MASTER
!$ IF (OMP_GET_NUM_THREADS() /= 0) THEN
!$    USE_OPENMP = .TRUE.
!$    OPENMP_AVAILABLE_THREADS = OMP_GET_NUM_THREADS()
!$ ENDIF
!$OMP END MASTER
!$OMP BARRIER
!$OMP END PARALLEL

END SUBROUTINE OPENMP_CHECK

END MODULE OPENMP

MODULE MISC_FUNCTIONS

IMPLICIT NONE

CONTAINS

SUBROUTINE SEARCH_CONTROLLER(NAME,CTRL_ID,DEVC_ID,DEVICE_INDEX,CONTROL_INDEX,INPUT_INDEX)

USE DEVICE_VARIABLES, ONLY: DEVICE,N_DEVC
USE CONTROL_VARIABLES, ONLY: CONTROL,N_CTRL
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
CHARACTER(100) :: MESSAGE
CHARACTER(*), INTENT (IN) :: NAME,CTRL_ID,DEVC_ID
INTEGER :: I, DEVICE_INDEX,CONTROL_INDEX
INTEGER , INTENT(IN) :: INPUT_INDEX
 
! There cannot be both a device and controller for any given entity

IF (DEVC_ID /= 'null' .AND. CTRL_ID /='null') THEN
   WRITE(MESSAGE,'(A,A,1X,I3,A)')  'ERROR: ',TRIM(NAME),INPUT_INDEX,' has both a device (DEVC) and a control (CTRL) specified'
   CALL SHUTDOWN(MESSAGE)
ENDIF

! Search for device

IF (DEVC_ID /= 'null') THEN
   DO I=1,N_DEVC
      IF (DEVICE(I)%ID==DEVC_ID) THEN
         DEVICE_INDEX = I
         RETURN
      ENDIF
   ENDDO
   WRITE(MESSAGE,'(A,A,A)')  'ERROR: DEVICE ',TRIM(DEVC_ID),' does not exist'
   CALL SHUTDOWN(MESSAGE)
ENDIF

! Search for controller

IF (CTRL_ID /= 'null') THEN
   DO I=1,N_CTRL
      IF (CONTROL(I)%ID==CTRL_ID) THEN
         CONTROL_INDEX = I
         RETURN
      ENDIF
   ENDDO
   WRITE(MESSAGE,'(A,A,A)')  'ERROR: CONTROL ',TRIM(CTRL_ID),' does not exist'
   CALL SHUTDOWN(MESSAGE)
ENDIF

END SUBROUTINE SEARCH_CONTROLLER


END MODULE MISC_FUNCTIONS

